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Hydrodynamic flow occurs in an electron liquid when the mean free path for electron-electron 
collisions is the shortest length scale in the problem. In this regime, transport is described by the 
Navier-Stokes equation, which contains two fundamental parameters, the bulk and shear viscosi¬ 
ties. In this Article we present extensive results for these transport coefficients in the case of the 
two-dimensional massless Dirac fermion liquid in a doped graphene sheet. Our approach relies on 
microscopic calculations of the viscosities up to second order in the strength of electron-electron 
interactions and in the high-frequency limit, where perturbation theory is applicable. We then use 
simple interpolation formulae that allow to reach the low-frequency hydrodynamic regime where 
perturbation theory is no longer directly applicable. The key ingredient for the interpolation formu¬ 
lae is the “viscosity transport time” r v , which we calculate in this Article. The transverse nature of 
the excitations contributing to r v leads to the suppression of scattering events with small momentum 
transfer, which are inherently longitudinal. Therefore, contrary to the quasiparticle lifetime, which 
goes as —1/[T 2 1ii(T/Tf)], in the low temperature limit we find r v ~ l/T 2 . 

PACS numbers: 73.20.Mf,71.45.Gm,78.67. Wj 


I. INTRODUCTION 

Hydrodynamics 1-3 is a powerful non-perturbative the¬ 
ory to deal with transport properties of strongly inter¬ 
acting many-particle systems. 

In the solid state, interactions need to be sufficiently 
strong to ensure that the mean free path £ ee = vpT ee 
for electron-electron (e-e) collisions is the shortest length 
scale in the problem, i.e. £ ee <C £ p ,L,vf/lo. Here, uf is 
the Fermi velocity, r ee is the quasiparticle lifetime due 
to electron-electron collisions 4 , l p is the mean free path 
for momentum-non-conserving collisions, L is the sample 
size, and co is the frequency of the external perturba¬ 
tion. In Fig. 1 we show the result of microscopic calcu¬ 
lations of the e-e mean free path for the two-dimensional 
(2D) massless Dirac fermion (MDF) liquid in a doped 
graphene sheet 1-9 embedded between two semi-infinite 
uniform and isotropic media with dielectric constants ei 
and 62 - The two-dimensional Fourier transform of the 
e-e interaction in this case is v q = 2 ire 2 / (eq), where 
e = (ei + 62)/2. Technical details on these many-body 
diagrammatic perturbation theory calculations can be 
found, e.g., in Refs. 5 and 6. We clearly see that, for 
sufficiently large temperatures, there is a wide range of 
carrier concentrations in which I ee becomes much shorter 
than the typical device size (L ~ 10 /im). 

The regime defined by the above inequalities is named 
in what follows “hydrodynamic”, “low-frequency” or 
“collisional”. In the range of temperatures and carrier 
densities in which this regime is attained, e-e interactions 
drive the system towards a local quasi-equilibrium state 
characterized by slowly-varying time-dependent density 
n{r,t) and drift velocity v(r,t), which obey the conti¬ 


nuity and Navier-Stokes equations. The latter are con¬ 
trolled by two transport coefficients, the shear viscos¬ 
ity, rju^o, which describes the friction between adja¬ 
cent layers of fluid moving with different velocities, and 
the bulk viscosity, Cw-> 0 j which describes the dissipation 
arising in the liquid when it undergoes a homogeneous 
compression-like deformation 1 . 

When, on the contrary, the frequency u> of the external 
perturbation is much larger than the quasiparticle colli¬ 
sion rate (i.e. WT ee 1)- but still much smaller than the 
characteristic free-particle frequencies epitomized by the 
Fermi energy—e-e interactions fail to drive the system 
towards local quasi-equilibrium 4 . Nonetheless, it is still 
possible to describe the system by hydrodynamic equa¬ 
tions of motion 10 , provided that the low-frequency bulk 
and shear viscosities are replaced by their high-frequency 
counterparts 4 (Coo and 7700, respectively) and that a finite 
value of the shear modulus (Sac) is allowed. This regime 
is named in what follows “high-frequency” or “collision¬ 
less” . We emphasize that it is a “high-frequency” regime 
only on the scale of e-e collisions, but not at all on the 
scale of the Fermi energy. 

In this Article we calculate the frequency-dependent 
viscosities and rj u for the 2D MDF liquid in a doped 
graphene sheet 8 . Doping, which creates a Fermi liq¬ 
uid of electrons or holes in the upper or lower Dirac 
band, is of essence here: we note that the zero-frequency 
shear viscosity of thermally excited electron-hole pairs in 
an undoped graphene sheet was previously calculated in 
Ref. 11. The low-frequency bulk and shear viscosities of 
an ordinary three-dimensional (3D) parabolic-band elec¬ 
tron gas in the Fermi liquid regime were calculated long 
ago by Abrikosov and Khalatnikov 12 . They found that 
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77w->o = 5oo7" v , where r v is a “viscosity transport time” 
of the order of (but not identical to) r ee (here all quan¬ 
tities refer to a 3D system). The Abrikosov-Klralatnikov 
calculation was extended to the opposite regime of high 
frequency in Ref. 10. The main difficulty in connect¬ 
ing these two regimes lies in the fact that e-e interac¬ 
tions play very different roles in the two cases. In the 
low-frequency regime their main effect is to cut off an 
otherwise diverging shear viscosity: the corrected shear 
viscosity is proportional to r v , which is non-perturbative 
in the strength of e-e interactions. Conversely, in the 
high-frequency regime e-e interactions generate a non¬ 
zero value for the otherwise vanishing shear viscosity: 
this finite value can and has been calculated perturba- 
tively. 

Our theoretical approach is based on ideas first pre¬ 
sented in Ref. 10. We combine the perturbative informa¬ 
tion contained in rjoo and £00 with the calculation of r v 
to generate non-perturbative interpolation formulas for 
riuj and The latter are approximately valid at all fre¬ 
quencies and consistently include many-body self-energy 
and vertex corrections. The final formulas are 10 : 


C , (^v ) 2 
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where, in addition to the above-mentioned quantities r 
and Coo, we also see the high-frequency shear modulus 
Soo, which is to the shear viscosity what the “Drude 
weight” is to the conductivity. We note that the shear 
modulus is renormalized by e-e interactions: however, 
this effect is relatively small, and it is thus qualitatively 
correct to approximate by its non-interacting value, 
which in the case of graphene is = nep/A, where £p 
is the Fermi energy. 

From a mathematical point of view, the viscosities ap¬ 
pear as coefficients in the expansion of the stress tensor 
t„„ to first order in the spatial derivatives of oscillating 
velocity fields, v^ = \ (d^v^ + d v v^). This can also be 
viewed 13 as the out-of-phase component of the response 
of to an oscillating metric field g MI ,. For a 2D isotropic 
fluid the expansion has the form 


V'M = (C w - r? w )[V • v(u})]5^ + 2q u v^ v (u) . (2) 


In a parabolic-band electron gas 4 , and also in graphene 
in the Fermi liquid regime 8 , the response of the stress 
tensor to the metric field is connected by equations of 
motion to the non-local response of the current to a vec¬ 
tor potential, i.e. the coefficient of q 2 in the expansion 
of the non-local conductivity for small wave vectors q. 
This implies that the high-frequency viscosities can be 
extracted from the damping rate of plasmons, the high- 
frequency collective excitations of an electron liquid 4 ' 14 . 
On the other hand, no standard protocol exists at present 
to measure the low-frequency viscosities of electrons in a 
solid-state host. Tomadin et al. 15 proposed a Corbino 


disk viscometer, which allows a determination of the hy¬ 
drodynamic shear viscosity rj u ^.0 from the dc potential 
difference that arises between the inner and the outer 
edge of the disk in response to an oscillating magnetic 
flux. More recently, it has been shown 16,1 ' that 77 w _>o 
can also be extracted from purely-dc non-local transport 
measurements in ultra-clean multi-terminal Hall bar de¬ 
vices. 

Our Article is organized as follows. In Sect. II we in¬ 
troduce the tight-binding model of graphene, which we 
use 18-20 to avoid broken gauge invariance due to the pres¬ 
ence of a rigid ultraviolet cut-off in the MDF low-energy 
theory. The MDF limit is indeed taken only at the very 
end of the calculation. In Sect. Ill we derive the relativis¬ 
tic counterpart of the Navier-Stokes equation 1 , which de¬ 
scribes the long-wavelength dynamics of quasiparticles in 
graphene in both the low- and high-frequency regimes. 
In this Section we also show the connection between the 
macroscopic bulk and shear viscosities and the longitu¬ 
dinal and transverse current-current response functions, 
which can be microscopically calculated from the usual 
Kubo formula 4 . In Sect. IV we use a kinetic equation ap¬ 
proach and the relaxation-time approximation to deter¬ 
mine the interpolation formulas for the viscosities and the 
elastic moduli of graphene in terms of the high-frequency 
viscosity 77 ^ and a yet undetermined transport time r. 
These quantities are calculated in the remainder of the 
Article. Since the longitudinal current-current response 
function of the 2D MDF liquid in a doped graphene sheet 
was calculated in Ref. 18, in Sect. IV A we focus on its 
transverse counterpart. We calculate it at the lowest non¬ 
vanishing order in the strength of e-e interactions, which 
is quantified by the graphene’s fine structure constant 8 
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Our results for the high-frequencies bulk and shear vis¬ 
cosities are reported in Sect. IV A. We prove that the 
high-frequency bulk viscosity vanishes, while 77 ^ is finite. 
Sect. IV B is devoted to the calculation of the viscosity 
transport time r v . The approach we adopt is very similar 
to that used in Refs. 21 and 22, where the e-e contribu¬ 
tions to the charge, spin, and thermal conductivities of 
the 2D MDF liquid in a doped graphene sheet were cal¬ 
culated. Therefore, only the main steps of the calculation 
are surveyed. We refer the reader interested in more de¬ 
tails to Ref. 22. Finally, In Sect. V we show our result 
for the shear viscosity at finite frequency 77 ^. Making 
use of Eq. (1), we provide numerical results for the shear 
viscosity of the 2D MDF liquid in a doped graphene at 
all frequencies. Appendix B presents a self-contained de¬ 
scription of the generalized relaxation time approxima¬ 
tion, leading to the formulas of Eq. (1). Appendix D 
contains several technical details of the calculation. In 
this Article we set, except when explicitly stated other¬ 
wise, h = 1 . 
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FIG. 1. (Color online) Panel a) The e-e mean free path 
fee (in /im) in a 2D MDF liquid is plotted as a function 
of temperature T (in K). Different curves refer to different 
values of the excess carrier density, n = 0.5 x 10 12 cm -2 , 
n = 1.0 x 10 12 cm -2 , and n = 2.0 x 10 12 cm -2 . Panel b) 
The e-e mean free path l ee (in /rm) in a 2D MDF liquid is 
plotted as a function of the excess carrier density n (in units 
of 10 12 cm -2 ), for three values of T, i.e. T = 100, 200, and 
300 K. In both panels the e-e dimensionless coupling constant 
has been set to a ee = 0.5. 

II. MODEL AND BASIC DEFINITIONS 

Following Refs. 18-20, we describe 7r-electrons in 
graphene by a one-orbital tight-binding (TB) model 1 . To 
keep the model as simple as possible, we set to zero all the 
hopping parameters but the nearest-neighbor one. The 
low-energy MDF limit will be taken only at the very end 
of the calculation, after carrying out all the necessary 
algebraic manipulations. By following this procedure we 
avoid problems associated with the introduction of a rigid 
ultraviolet cut-off, which breaks the gauge invariance 23 
and is responsible for the appearance of anomalous com¬ 
mutators 23 ’ 24 . 

The non-interacting Hamiltonian is 

Tio= Y ' <r<xp)^k,t> , ( 4 ) 

fe6BZ,a:,/3 

where the operator ^ a (‘tpk,a) creates (annihilates) an 
electron with Bloch momentum k, which belongs to the 
sublattice 7 a = A, B. The vector fk is defined as' 

3 

fk = (He [e~ lk ' 5 '} , -3m [e ~ lk - s '}) . (5) 

i -1 

Here t ~ 2.8 eV is the nearest-neighbor tunneling am¬ 
plitude, while Si (i = 1,...,3) are the vectors which 


connect an atom to its three nearest neighbors, i.e. 
<$i = aV3x/2 + ay/2, S 2 = —ax/ix/2 + ay/2, and 
£3 = —ay. Here a ~ 1.42 A is the Carbon-Carbon 
distance in graphene. The sum over k in Eq. (4) is re¬ 
stricted to the first Brillouin zone (BZ) and the Pauli 
matrices cr l a p (i = x,y,z ) operate on the sublattice de¬ 
gree of freedom. 

The TB problem posed by the Hamiltonian (4) can be 
easily solved analytically'. One finds the following eigen¬ 
values £k,x = A|/fe|, with A = ±. These two bands touch 
at two inequivalent points (K and K') in the hexago¬ 
nal BZ. The low-energy MDF model is obtained from 
Eq. (4) by taking the limit a —> 0, while keeping the 
product ta constant. In this limit fx+k —> 'C’f k, where 
= 3fa/2 ~ 10 6 m/s is the density-independent Fermi 
velocity. 

Introducing the field operator cj, A (ck,x) as the cre¬ 
ation (annihilation) operator in the eigenstate represen¬ 
tation, Eq. (4) can be rewritten as 

iio = A Cfc,A • (6) 

k,\ 


In the same representation the Hamiltonian describing 
e-e interactions reads 4 


ii 


ee 
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where the density operator is 
n q = P AA'( fe - + 9/2)4- g /2,A fi fc+q/2,A' » 

k, A,A' 

( 8 ) 


and v q is the 2D discrete Fourier transform of the real- 
space Coulomb interaction, which is a periodic func¬ 
tion of the reciprocal-lattice vectors, and reduces to 
~ 2ire 2 /(eq) in the limit of q —> 0. Finally, in Eq. (8) 
we have introduced the “density vertex” 

e i(e k -e h ,)/2 , \\’ e -i( 8 h -e kl )/2 
V X \' ( k , k') = - - - (9) 

with 9 k = Arg [f k , x + ifk, y }- Here {fk,i,i = x, y} denotes 
the Cartesian component of the vector fk- In the low- 
energy MDF limit, 9x+k —t fk, where fk is the angle 
between k and the x axis. 

Note that in writing Eq. (7) we have neglected the one- 
body operator proportional to the total number of parti¬ 
cles, which avoids self-interactions 4 , since it has no effect 
on the calculations we will carry out below. The viscosi¬ 
ties are indeed determined (at the lowest non-vanishing 
order in the strength of e-e interactions) by two-particle 
excitations only, which are generated by two-body opera¬ 
tors. 

We also introduce the current operator 

Jq,a =Y Y ~ <?/ 2 ’ k + 9 / 2 ) 

k,/3 A,A' ^ 

X ^1-<j/2,aD s +9/ 2 A ’ (I®) 
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where the “pseudospin-density” vertices are 


, . V e i(S k +9 k ,)/2 , \ e -i(8 k +8 h ,)/2 

s$,(k,k') = — -- ( 11 ) 


and 


S™(k,k') = 


y e i{e k +e k ,)/2 _ Xe -i(e k +e k ,)/2 
2 i 


■ ( 12 ) 


For the sake of definiteness, we assume the system to 
be ro-doped, with an excess electron density n. Results 
for a p-doped system can be easily obtained by appealing 
to the particle-hole symmetry of the MDF model defined 
by Eqs. (6)-(7). The Fermi wave vector is defined as 
ftp = y47 mjN{, while £f = vf^f is the Fermi energy, 
and Nf = 4 is the number of fermion flavors in graphene. 
For future purposes we also define the matrix element of 
the z component of the pseudospin between the states 
labeled by fe, A and k',X': 




gi(0k Q k r)/2 — c 
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III. THE BULK AND SHEAR VISCOSITIES — 
GENERAL THEORY 

The bulk and shear viscosity are usually introduced 1 
as phenomenological coefficients to describe the long- 
wavelength motion of a viscous fluid close to a quasi¬ 
equilibrium situation. The interactions between the el¬ 
ementary constituent of the fluid, although extremely 
complicated at the microscopic level, admit a rather sim¬ 
ple description in terms of macroscopic coefficients. Their 
space and time average is in fact responsible for the fric¬ 
tion between (macroscopic) fluid elements having differ¬ 
ent values of the momentum. The two viscosities then 
describe the forces between fluid elements that undergo 
either a shear or a compression-like long-wavelength de¬ 
formation. 

It is therefore possible, starting from a fluid-like de¬ 
scription of the 2D MDF liquid in graphene, to derive the 
macroscopic response of the current to an external vec¬ 
tor potential, in the linear response regime and in terms 
of hydrodynamic coefficients. Equating the coefficient of 
proportionality between the current and the vector po¬ 
tential to the microscopic current-current linear response 
functions of the system, one obtains a microscopic defi¬ 
nition of the bulk and shear viscosities. It turns out that 
these can be calculated from the coefficients of the expan¬ 
sion to order q 2 /uj 2 of the current-current linear response 
functions (in the limit q —> 0). This approach was used 
in Ref. 10 for the case of a parabolic-band ( i.e. Galilean 
invariant) electron gas. 

A more fundamental approach 25-27 relies on the fact 
that it is possible to microscopically define the stress ten¬ 
sor operator starting from the equation of mo¬ 

tion of the momentum density, and to calculate the de¬ 
viation of its average from the equilibrium value due to 


an applied strain. In the linear regime and at low fre¬ 
quencies, such variation is proportional to the applied 
strain. The coefficient of proportionality is the “tensor 
of elasticity”, a rank-4 tensor whose imaginary part in a 
rotationally and time-reversal invariant system at q = 0 
can be characterized by two coefficients, the complex bulk 
and shear moduli. Note that in the linear regime the ten¬ 
sor of elasticity is equivalent to the stress-stress response 
function. Therefore, the knowledge of the latter response 
function constitutes a viable route to the calculation of 
the viscosity coefficients 25-27 . 

The connection between the two approaches is rather 
trivial in Galilean invariant systems. Indeed in such sys¬ 
tems the current density is proportional to the momen¬ 
tum density, and therefore its time derivative is propor¬ 
tional to the divergence of the stress tensor. It is there¬ 
fore possible to derive an equation of motion that relates 
the current-current response functions to the stress-stress 
response. From that, it is thus evident that the viscos¬ 
ity, which is proportional to the coefficient of the term 
of order q 2 /ai 2 in the expansion of the current-current 
response function, can also be calculated from the q = 0 
limit of the stress-stress response. 

In the 2D MDF liquid, however, the current and the 
momentum are not proportional to each other. The cur¬ 
rent is indeed an off-diagonal operator (in pseudospin 
space), which represents the hopping between the two 
inequivalent Carbon atoms in the unit cell. Conversely, 
the momentum operator is a diagonal one (as, e.g., the 
density operator). Therefore the two approaches could 
give, in principle, different results. Nevertheless, we find 
that in the Fermi liquid regime, i.e. when graphene is 
doped and the temperature T <C Ef/^’b is sufficiently 
low (&b is the Boltzmann constant), the two agree with 
each other. In the following subsections we show in de¬ 
tail how the two approaches reconcile in the Fermi liquid 
regime. Here, however, we briefly discuss the problem in 
general terms. 

It is quite easy to understand what goes wrong in the 
general case. While the second approach, i.e. the calcu¬ 
lation from the stress-stress response, is completely gen¬ 
eral, the first one relies on the fact that it is possible to 
write a Navier-Stokes equation for the velocity as in the 
classical non-relativistic case. This statement is highly 
non-trivial. Indeed, in graphene the velocity operator, 
being proportional to the unit vector of the momentum, 
is not a conserved quantity. Therefore, it is not possi¬ 
ble to write in general a Navier-Stokes equation for the 
macroscopic velocity in the same way as it is done in 
the Galilean invariant case. Since the quantities that are 
conserved by e-e interactions are the density, energy den¬ 
sity, and momentum, one should in principle consider the 
Navier-Stokes equations for these three. In passing we 
note that, in a system with a linear energy dispersion, 
the momentum coincides, up to a proportionality con¬ 
stant, with the non-interacting component of the energy 
current. 

However, as it was shown in Ref. 22, in the Fermi liq- 
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uid regime the velocity is essentially a conserved quantity. 
This statement is rationalized as follows. At low temper¬ 
ature and in a doped system, the states that contribute 
to transport are those in a narrow region (of size ~ kp,T) 
around the Fermi energy. Those states have momentum 
nearly identical to the Fermi momentum kp, and the ve¬ 
locity operator evaluated on them is simply proportional 
to the momentum (via the constant kp ). Since the same 
states give the dominant contribution to the total mo¬ 
mentum of the system, and the latter is conserved by e-e 
interactions, also the velocity is a conserved quantity 28 . 
This paves the way to a hydrodynamic description of the 
velocity, with a Navier-Stokes equation 29 that contains 
the same viscosities as the momentum one. It is indeed 
clear that the two equations must be proportional to each 
other via the cyclotron mass m c = kp/vp. 


A. The hydrodynamic approach 

The long-wavelength dynamics of the 2D MDF liquid 
in a doped graphene sheet is described by the relativis¬ 
tic Navier-Stokes equation in 2+1 dimensions 1 . In or¬ 
der to avoid confusion, we remark that the use of a rel¬ 
ativistic equation is dictated by the low-energy linear- 
in-momentum band dispersion of quasiparticles, which, 
however, move with a Fermi velocity vf that is much 
smaller than the speed of light c. For this reason we 
neglect retardation effects. 

The Navier-Stokes equation for MDFs reads 29 

wu^d p u v - (8£ - u^u^dfj.P + (6% - u v u^)d p T^ = F u , 

( 14 ) 

where the space- and time-dependences of the various 
quantities are suppressed for brevity and the sum over 
repeated indices is understood. Here a M is a contravari- 
ant vector, while a p = rj pv a v is its covariant counterpart, 
and = diag(l, —1, —1, —1) is the metric tensor of flat 
spacetime. 6 v is the four-dimensional Kronecker delta. 
In Eq. (14) w = w(r,t) is the enthalpy density (which is 
equivalent to the Drude weight), P = P(r , t) is the pres¬ 
sure, i+ = 7 (i>)(l, v/vp) is a three-component velocity 
vector, v = v(r,t) is the 2D fluid-element velocity, and 
7 (v) = (1 — v 2 /vp)- 1 / 2 . We assume that all the ther¬ 
modynamic quantities in Eq. (14) have been regularized 
by subtracting the (divergent) non-interacting contribu¬ 
tion due to the infinite sea of states in valence band. 
This implies that in the non-interacting limit they de¬ 
pend only on the excess carrier density in the conduction 
band. However, corrections due to e-e interactions be¬ 
tween states in conduction and in valence band are fully 
taken into account. 

Moreover, we assume that states in valence band do 
not contribute to transport, which is certainly a good ap¬ 
proximation in the limit in which the wave vector and fre¬ 
quency of the external perturbation are small. The fluid 
velocity can thus be expressed, in the “non-relativistic” 
limit |u| <C vf, in terms of the macroscopic current 


j(r,t) carried by quasiparticles in conduction band as 
v(r,t) = j(r,t)/n(r,t). Here n(r,t) is the position- and 
time-dependent excess number density. In Eq. (14) we 
also defined the stress tensor r= r M ^(r, t), whose most 
general form, based on symmetry arguments, is reported 
in Eq. (2) 1 . Finally, in the Navier-Stokes Eq. (14) we in¬ 
troduced the driving force 4 F v = F„{r,t). As explained 
above, we do not introduce any retardation effect in the 
driving term since the Fermi velocity is much smaller 
than the speed of light. 

Linearizing Eq. (14), taking the “non-relativistic” limit 
|u| -C tip [which implies q(v) ~ 1], and considering only 
the spatial components we get 1,11 ’ 29 

~2dtVi + Vi + -kd t ) P-V A r ) 0 [V^ + VjU,- 

V F \ V F / 1 

- (V • v)5ij] + Co(V • w)<5jj| = n^dt.Ai , (15) 

where the Latin indices i,j = x,y denote the Cartesian 
components of the vectors. Since we are interested in 
deriving the current response of the fluid to an external 
perturbation, in Eq. (15) we expressed 4 the driving field 
as the time derivative of a vector potential A = A(r,t)- 
We assume the density to be close to its equilibrium 
value (which only in this Section we denote as n eq ), 
i.e. n(r,f) = n eq +5n(r,t) with \Sn(r, t)/n eq \ <C 1. More¬ 
over, in a linear-response fashion, also the current density 
j(r , t) and the vector potential A(r , t) are assumed to be 
small (note that the current is zero at equilibrium). 

Using well-known thermodynamic relations we 
rewrite 4 

dP B 

d t P = -p—d t n =- d t n , (16) 

on n eq 

and 

dP B 

VP=—Vn= —Vn, (17) 

On n eq 

where B = B(n) is the bulk modulus 4 of the 2D MDF 
liquid regularized, as explained above, by subtracting 
the non-interacting contribution of the valence band. 
In Eqs. (16) and (17) we suppressed the space- and 
time-dependences of the various quantities for brevity. 
These equations, together with the continuity equation 
d t n(r,t) = —V • j{r,t), allow us to rewrite the Fourier 
transform of Eq. (15) as 

u) 13 & 

-2^ji--{Q-j)Qi+iVoQ 2 3i+Ko(Q-j)Qi = n 2 -ujA z . (18) 

To obtain this equation we discarded a term proportional 
to the product 5n(r, t)A(r, t), which gives a contribution 
beyond the linear regime. Moreover we used that in this 
regime d t v(r,t ) = d t \j{r,t)/n(r,t)] ~ d t j(r,t)/n e q . A 
similar approximation was used for terms which contain 
spatial derivatives of v(r,t). 

We now observe 4 that Eq. (18) is suitable to describe 
the long-wavelength dynamics of 2D MDFs in the col- 
lisional uiT ee <C 1 limit. In this regime 4 the system is 
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indeed expected to behave as a liquid, with a finite shear 
viscosity and a vanishing shear modulus. However, in the 
opposite, collisionless, limit (wr ee 1) a solid-like behav¬ 
ior is expected to emerge, characterized by a small shear 
viscosity and a finite shear modulus 4 . The system may 
thus be described by the equations of the elasticity the¬ 
ory 30 . A unified description is achieved by introducing 4 
the complex bulk and shear moduli, B u = B u — icvCu and 
Sui = S u — iojrju, respectively, and by rewriting Eq. (18) 
as 

—ujji - (q ■ 3 )q i - q j z = n -uAi . (19) 

Up OJ OJ H C 

Note that, from the above discussion we expect 
3?e(5 w ) —> 0 in the low-frequency regime. Separating 
the longitudinal and transverse components of Eq. (19), 
we finally get 


the two viscosities can be calculated perturbatively, the 
first non-vanishing contribution being of second order in 
the dimensionless strength of e-e interactions, i.e. a ee in 
Eq. (3). It is therefore legitimate, up to second order, to 
neglect interaction corrections to w. 

By inverting Eqs. (23) and (24) we can express the 
high-frequency bulk and shear viscosities of a 2D MDF 
liquid in terms of the current-current response functions. 
To second order in the strength of e-e interactions these 
hydrodynamic coefficients read 4 

T ]oo = - lim lim , 

u — >-0 q — >-0 q z 

2 

(jjrn 

Coo = - lim lim —^ {Sm[x L (g,u)] - 3 w[xt(<7,w)]} . 

Ul —>0 q —>0 q z 

(25) 


jh(q^) 


2/ 2 - ~—-A L (q,w) 

wbJ 2 /v p - {B u + S u )q 2 C 

XL(g,w)^A L (g,w) , (20) 


and 


3 t (q,u>) = 


2 2 


J 2 /v 2 - S u q* c 


A T (q,u) 


= XT (q,bj)-A T (q,ui) . 


( 21 ) 


The transverse component of any 2D vector v is defined 
as dt = v — q{q ■ v), while its longitudinal component 
is ul = q ■ v. Eqs. (20) and (21) provide a macroscopic 
definition of the longitudinal [xl(< 3S w )] and transverse 
Lyt(<Z,w)] components of the current-current linear re¬ 
sponse function. Indeed, in a rotationally-invariant sys¬ 
tem 


/ \ / \ QtQk 

XuihVhu) = XL (q,u>)-^- 

+ Xr{q,u) ($ik ~ i (22) 

where t, k = x,y are Cartesian indices. It is easy to show 
that, if q = qx, Xl (q,v) = Xj x jM^) and XT (g,w) = 
Xj y j y (q,bj). These relations will be used in what follows. 
Expanding Eqs. (20) and (21) to order q 2 /u> 2 we get 

( 2 \ 2 2 

^|y^ eq j ^jiCu+Vu), (23) 


and 


3m[x T (g,w)] 




(24) 


where we have approximated the enthalpy density by 
its non-interacting value = n eq £F- This is justi¬ 
fied because the relations (23)-(24) will be used to esti¬ 
mate the high-frequency bulk and shear viscosities, Coo 
and 77oo, respectively. In the high-frequency regime, 


Note that in Eq. (25) the limit u> —> 0 is taken after the 
limit wr ee ^ 1. 

The current-current response functions on the right- 
hand side of Eq. (25) have a rigorous microscopic defini¬ 
tion 4 in terms of Kubo products, i.e. 31 

Xab{u) = ^((A;B)) u , (26) 


where S is the 2D electron system area, A and B are 
operators, and 

/»oo 

{{A-,B)) u = -i dte^+^dAq),^}) . (27) 

Jo 

Note that the average (...) in Eq. (27) is taken over the 
ground state of the interacting system. 

The longitudinal component xl (<?,w) of the current- 
current response function of 2D MDFs was calculated 
in Ref. 18. In this Article we evaluate the transverse 
current-current response function XT(q,oj) at second or¬ 
der in the strength of e-e interactions and in the limit 
vpq «w« 2s F - 

We note that, from Eq. (20), it is possible to derive an 
expression for the non-local charge conductivity a(q, w). 
Replacing w 2 —> bj[bj + i/r(q,bj)\ in its denominator to ac¬ 
count for non-momentum-conserving dissipative effects, 
and neglecting the real parts of B u and S u (which are 
negligible in the limit %g <C w) we get 


er(g,u;) 


Xl (g,uQ 
—iuj 


n/m c 

—itjJ + l/r(g, to) + v^q 2 


■ (28) 


In Eq. (28) we defined the kinematic viscosity = 
Vui/{n-m c ), and we used the fact that Cu> = 0 (as will 
be demonstrated in what follows). Since Eq. (28) is valid 
to all orders in the perturbative expansion in the strength 
of the e-e coupling constant, the cyclotron mass m c must 
be renormalized according to Landau theory of normal 
Fermi liquids 4 . 
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B. The microscopic approach 


In order to get a macroscopic Navier-Stokes equation 
for a given quantity, it is necessary for the latter to be 
conserved by interactions at the microscopic level. For 
example the momentum density operator p(x,t), whose 
Fourier transform is microscopically defined as 

Pq = a > (29) 

fcGBZ,Q: 


is a locally conserved quantity, and satisfies the following 
continuity equation 


d t pj{x,t) = - diTij(x,t) , (30) 


where i,j = x,y are Cartesian indices and the summa¬ 
tion over repeated indices is understood. Eq. (30) defines 
the symmetric stress tensor operator Tij{x. t). In the 
non-interacting limit, the Fourier transform of fij(x , t) is 
defined by 


do) 


?£ = <*£ 3 


t 

k—q, a 


kia° n o + hi<ri 


■ ( 31 ) 




Let us now consider a time dependent deformation 
u(r,t) of the electron coordinates {r l: i = 
which are thus shifted as r* —► r,; — u{ri,t). This de¬ 
formation should not be confused with a similar one that 
could be applied to the lattice. To linear order, the met¬ 
ric tensor of the deformed system is 

Qij — $ij ( D,Uj T djUi) 

= Sij - 2 Uij , (32) 

where Ujj is the strain tensor of the elasticity theory 30 . 
Note that g 13 = S 13 + 2 u 13 . Under this transformation 
the Hamiltonian becomes (to linear order in utj ) 


H' =H + 2 — 
6gv 

9 

= H+ TijU lj . 


(33) 


It can be shown 32 that the non-interacting part of the 
stress tensor defined by Eq. (33) is identical to Eq. (31). 
The variation of the Hamiltonian, TijU 13 , induces a varia¬ 
tion in the expectation value of the stress tensor operator. 
To linear order in u 1 ' 3 we get 


S(Tij)(q,uj) = Q. ljM {q,u})u kl {q,u}) , (34) 


which defines the tensor of elasticity Qij t kl(q,w). In a 
rotationally invariant system this tensor can be decom¬ 
posed as 


lim Qij,ki(q 

uj —>-0 


0 , oj) = B u 5ij6u 






( 35 ) 


where d = 2 is the dimensionality of the system. From 
the general theory of linear response 4 and the form of 
the perturbation in Eq. (33), it is clear that the elasticity 
tensor can be calculated from the stress-stress response 
function, i.e. 

Qij, kite ’ w ) = jg{{Tij{q);Tki(-q))) u ■ (36) 

In writing Ecp (36) we have neglected a “contact” term 
(analogous to the diamagnetic term of the current- 
current response function), which has been recently dis¬ 
cussed e.g. in Ref. 33. However, this term is purely real 
and thus does not contribute to the bulk and shear vis¬ 
cosities, which are defined as 


rju = - lim 


%m[Xxy,xy{q = 0,w)] 


Cw = - lim • 

uj —>-0 


^rn[xxx,xx(q = O.qQ] 

UJ 


2 {d- l)7sm[xxy,xy(q = o,w)]/d 


(37) 


C. Equivalence of the hydrodynamic and 
microscopic approaches in a Fermi liquid 

In a Galilean invariant system (in which the current 
operator is j q = p q /m c ) it is rather easy to see that the 
hydrodynamic and microscopic approach give identical 
results. Thanks to Eq. (30) and the following identity 

u;((A-B)) u = ((id t A-,B)) u + ([A,B)) 

= ((A--id t B)) u + ([A,B]) , (38) 

we get (q = qx) 

s m[xxy,xy(q = 0,w)] = lim -^-3m[x T (9,w)] , 

m 2 uj 2 

%rn[xxx,xx(q = 0,w)] = lim -^-3m[x L (g,w)] ■ 

(39) 

Here, we used twice Eq. (38) and the fact that the av¬ 
erage of the commutator on the right-hand-side of that 
equation is purely real. From Eqs. (37) and (39) one 
immediately recovers Eq. (25). Notice that, since in a 
Galilean invariant system m c is not renormalized by e-e 
interactions, Eq. (25) is valid to all orders in the strength 
of e-e interactions, and, therefore, at all frequencies (as 
long asw< £f)- 

The connection is much less straightforward in the case 
of the 2D MDF liquid in graphene. Indeed, in a system 
with a linear dispersion, the current operator is not pro¬ 
portional to the momentum-density operator. Therefore, 
Eq. (39) does not hold, in general. This fact notwith¬ 
standing, we have found that in the Fermi liquid regime 









Eq. (39) still holds in an approximate sense and, there¬ 
fore, to the level of accuracy we are interested in, the 
microscopic and hydrodynamic approaches coincide. The 
details of the argument are presented in Appendix A, but 
the basic idea is quite simple. First of all, it is clear that 
at sufficiently low frequency interband transitions are ir¬ 
relevant and can be disregarded. Then it is clear that, 
within a narrow band of energies around the Fermi level 
in the conduction band, there is no qualitative difference 
between the linear dispersion of MDFs and the linearized 
dispersion of ordinary massive fermions: the two disper¬ 
sions are indistinguishable if the ordinary Schrodinger 
fermions are assigned the cyclotron mass m c = kp/vF- 
We conclude that the equivalence of the microscopic and 
hydrodynamic approaches carries over to MDFs with the 
simple replacement of the effective mass by the cyclotron 
effective mass. 


IV. CALCULATION OF THE VISCOSITY IN 
THE RELAXATION TIME APPROXIMATION 


As we pointed out in the Introduction, e-e interactions 
enter the calculation of the viscosities quite differently in 
the high-frequency (collisionless) and low-frequency (col- 
lisional) regimes. In the high-frequency regime, the ef¬ 
fect of the interaction is perturbative, meaning that there 
would be no viscosity without interactions creating a cor¬ 
relation between the motions of adjacent parts of the liq¬ 
uid. In the low-frequency regime, e-e interactions are 
non-perturbative as their primary role is to establish a 
finite mean free path £ ee for electrons: this mean free 
path would be infinite in the absence of interactions. The 
problem is how to connect in a seamless way these two 
very different regimes of e-e scattering. The relaxation 
time approximation, summarized in Eqs. (1), offers a sim¬ 
ple and physically motivated way to achieve this connec¬ 
tion. The derivation of these formulas closely parallels 
the derivation given in Ref. 10 for Galilean invariant sys¬ 
tems. Only a minor adaptation is needed, namely the 
replacement of the bare electron mass m by the cyclotron 
effective mass m c , as discussed in the previous Section. 
We therefore refer the reader to Ref. 10, where a de¬ 
tailed derivation of Eqs. (1) is provided, and we focus in 
this Section on the calculation of the inputs for Eqs. (1), 
i.e. the high-frequency viscosities and the corresponding 
relaxation times. An alternative derivation of Eqs. (1) is 
presented in Appendix B. 


A. The high-frequency viscosities 

To evaluate the transverse current-current response 
function, we adopt the (Hamann-Overhauser 34 or 
Schrieffer-Wolff 35 ) canonical transformation approach 
outlined in Refs. 18-20. 

First of all, we introduce a unitary transformation gen¬ 
erated by a Hermitian operator F, 

n' =e lP (Ho + n ee )e- iP , (40) 

which cancels the e-e interaction from the transformed 
Hamiltonian, i. e. we require T~L' = "Ho to second order in 
the dimensionless strength a ee of e-e interactions. The 
transformation F = 1 + F\ + F 2 +... is determined order- 
by-order in perturbation theory. Here 1 is the identity 
and F n is the term of n-th order in the strength of e-e 
interactions. As shown in Refs. 18-20, to calculate the 
imaginary part of the current-current response function 
in the limit vpq w <C 2ep it is sufficient to determine 
F\, which satisfies the equality [H 0 , iF\\ = H ee . 

After carrying out the transformation F, the Kubo 
product in Eq. (27) is reduced to the evaluation of a 
non-interacting response function oc ((A 1 , B'))o tU . The 
subscript “0” here means that the average (...) is now 
performed over the ground state of the non-interacting 
system and that the time evolution is generated by 
Ho- However, the operators A! = e lF Ae~ lF and B' = 
e lF Be~ lF are now dressed by e-e interactions. The key 
idea is to realize that the calculation of c Am[xj a j fl (q,uj)] 
to second order in the strength of e-e interactions and in 
the limit zip q <C w -C 2ep requires only the knowledge 
of the transformed current-density operator j' to first 
order 18-20 , i.e. j\ q = j q +ji, q , where 

Jl.q = [iFi,j q ] • (41) 

Thus, to second order in the strength of e-e interac¬ 
tions, and for z;pg <C oj -C 2ep, we get the follow¬ 
ing exact-eigenstate (Lehmann) representation 4 of the 
current-current response function 

= -K^(Q\h,qArn)(rn\ji ,-^1°) 

m 

x 6(u> - U} m0 ) . 

(42) 

The calculation of F\ and j-\ , q is carried out in Appen¬ 
dices C and D. 

Here we quote the final formula for the two components 
(longitudinal and transverse) of the current-current re¬ 
sponse function, exact to second order in e-e interactions 
and in the large- N{ limit, which is 


9 


%m[xe{q,uj)\ = - 


E 

a ,( 3 = x,y 1 


d-q' 2 
(2tt y Vq ' 


dui' 


{r^ (q, q’)T { p {-q, -g')9fm[x^(«'> a/)]3m[x^E (?',<*>- a/)] 




( 43 ) 


In this equation i = L,T and xin (?,<*>), Xj°^ (?><*>), 

and Xnj^iq^) are the non-interacting density-density, 
current-current, and density-current response functions 
of a 2D gas of MDFs. The quantities {Ta\q,q');a = 
x, y\t = L, T} are defined as 


ri T W)=^ 


<lWy q’ a 


q' 2 k F 
Qx^&iy d - qy^OL,y 


- 1- 


4fc 2 


4v F k 3 F 6a ’ x ’ (44) 


and 


ri L) (<?,<?')= 



(45) 


We stress that the imaginary parts of the three linear- 
response functions Xnn(q,w), Xj°j> (?,<*>), and XnJ a (q,v) 
do not depend on any ultraviolet cut-off in the low-energy 
MDF limit. Moreover, since in the limit of w —> 0 the 
integral over q' is naturally restricted to 0 < q' < 2fc F , 
no ultraviolet regularization is needed in Eq. (43). The 
only pathology of the integral in Eq. (43) appears in the 
infrared q' —► 0 limit, due to the l/q' singularity of the 
Coulomb potential v q ’. This problem is cured by screen¬ 
ing, as we will further discuss below. 

We observe that, contrary to what happens in a 
parabolic-band electron gas, the matrix elements of 
Eqs. (44) and (45) do not vanish in the limit q —> 0. 
The terms that remain finite are due to broken Galilean 
invariance 23 , i.e. due to the presence of the valence band 
and, as noted in Ref. 18, are responsible for a finite 
optical conductivity in the single-particle optical gap 
u) < 2ep, which scales as ~ oj 2 . Being a conductivity, 
it is conceptually wrong to include it in the viscosities, 
and therefore in what follows we neglect the terms of 
Eqs. (44) and (45) that do not vanish when q —»■ 0, i.e. 
the last terms in Eqs (44)- (45). By insisting in retain¬ 
ing these terms, we would (wrongly) get two diverging 
viscosities. 

We stress that such a finite optical conductivity is 
an effect beyond the Fermi liquid theory, and that it 
is present only in the high-frequency limit. Note how¬ 
ever that, contrary to the derivation of the low-frequency 
limit, the calculation of the high-frequency viscosities 
does not require any Fermi-liquid assumption. There¬ 
fore, the general discussion of Sect. Ill does not need to 
be amended to take care of the peculiarities of the high- 
frequency calculation. 


Eq. (43) is the main result of this Section. Note that 
in the large-W limit the second-order expression for the 
imaginary part of the current-current response function 
in Eq. (43) has the appealing form of a convolution of 
two single-particle spectra 36 . The physical interpreta¬ 
tion of Eq. (43) is the following. At long wavelengths 
and to the lowest non-vanishing order of perturbation 
theory, the spectrum of the current-current correlator is 
dominated by the emission of two correlated electron- 
hole pairs. Each of the Kubo products on the right-hand 
side of Eq. (43) describes the rate of generation of a sin¬ 
gle electron-hole pair. The spectral weight associated 
with the excitation of two particle-hole pairs with oppo¬ 
site momenta and given total energy oj is proportional to 
their convolution. 

We now sketch how to make analytical progress in 
the evaluation of the current-current response function 
as from Eq. (43). 

The integrals in Eq. (43) can be carried out analyti¬ 
cally with the help of known formulas for the response 
functions 14 . We first observe that in the low-energy MDF 
limit the system is translationally and rotationally invari¬ 
ant. The current-current response function X'E/j^’ w ) 
is a rank-2 tensor that can be therefore decomposed ac¬ 
cording to Eq. (22). In the same way, the density-current 
response function can be seen to be equal to 

3m Xj°l(<? , i w ) = , (46) 


where Xj^ n Wi 1 ^) non-interacting longitudinal- 

current-density response function. Eqs. (22) and (46) can 
be used to perform analytically the angular integration 
in Eq. (43), which reads 



where 


a T = aL 


vW (g'VfcF ~ 4) 2 

w 4 32 


(48) 
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br. = 


v 2 q 2 3 q' A - 20k 2 q' 2 + 44A; 4 


U) 

vW ( q' 2 /k 2 - 2) 2 


32/cp 


&t = 


and, finally, 


or 


4 q 2 3</ 4 - 20ft 2 </ 2 + 44fcp 


32 


CL = 


CT — 


32fc 4 


CO 

vW W 2 /k 2 - 2) 2 

32 


(49) 


(50) 


To cure infrared divergences associated with the long- 
range tail of the e-e interaction v q i, we have eval¬ 
uated the integral over q' in Eq. (47) by replacing 
v q ' with a statically-screened Thomas-Fernri interaction, 
i.e. = 27 t e 2 /[e(q' + Qtf)]- Here <7 tf = N-pa ee kY is 

the Thomas-Fermi screening wave vector 8 . 

In the limit oj — > 0 we can expand the imaginary part 
of each response function Xrm, Xnj^i Xl°’ > ’ and Xt^ 
a power series of w', w — w' and retain only the lead¬ 
ing order of this expansion. The leading contribution to 
^mlxeiq, w)] in Eq. (43) in powers of w in the limit ui —> 0 
can be extracted from the following asymptotic formulae: 


hm n Sm[x[°i(g>')] = -N f 


u ' —>-0 


y^fcg-q' 2 u>' 

2 q' 2nvp 


(51) 


and 


co 


Jim SJmkfVV)] = ^T F ■ < 52 > 


The imaginary parts of density-current and lon¬ 
gitudinal current-current response functions scale 
with higher powers of a/. Indeed they can be 
derived from Eq. (51) according to the formu¬ 
lae 4 = a/Sm[xnn(gW)]/<?' and 

Sto[xl°VV)] = w ,2 3m[xin(g , ,w , )]/9 ,2 ■ 

Since the last two lines of Eq. (47) give a subleading 
contribution, in the limit v^q <C w 2ep, they can be 
disregarded. Plugging Eq. (47) back into Eq. (25) we 
finally get 

Voo — Ti*A.Nf (rr ee ) , 

Coo = 0 , (53) 

where ^Ar f (a e e) = 2 N { a 2 e f(Nfa ee ), with 


/O) 


15a: 3 - 15a: 2 - 52a: + 42 


(5a: 4 


288tt 

24a: 2 + 16) 
967T 


arccoth(l + x) . 


(54) 


Note that ^4Ar f (a ee ) is defined with an extra factor of 
2/lVf with respect to Ref. 18, due to the factor m 2 in 
Eq. (25), which has been absorbed into the definition 
of these functions. We underline that the dependence 
on the strength of e-e interactions beyond cn 2 e encoded 



FIG. 2. (Color online) The high-frequency shear viscosity 
r/oa (in units of the excess carrier density ?i) of the 2D MDF 
liquid in a doped graphene sheet is plotted as a function of the 
dimensionless parameter a ee , which defines the strength of e-e 
interactions. Note that the vertical axis must be multiplied 
by a factor 10 -2 . 


in f{x) is due to the Thomas-Fermi screened interaction 
introduced above used to regularize the infrared behavior 
of the integrand in Eq. (47). As expected 4 , the bulk 
viscosity is zero. 

In Fig. 2 we show the high-frequency shear viscosity, 
as defined in Eq. (53), in units of the carrier density n, 
and as a function of the dimensionless parameter a ee - In 
the limit of a ee —> 0 we find that 

N{ 

Voo -n—a 2 e ln(a ee ) . (55) 

07T 

The extra logarithmic dependence on the coupling con¬ 
stant of Eq. (55) is due to the Thomas-Fermi screening 
used in the calculation. For a ee 1 the Thomas-Fermi 
screened interaction becomes independent of a e e, and the 
high-frequency shear viscosity tends to —> n/(9nN{). 


B. The viscosity transport time 

We now turn to the evaluation of the viscosity trans¬ 
port time r v , which enters Eq. (1). This quantity can be 
estimated from a diagrammatic calculation of the low- 
frequency viscosity. Contrary to its high-frequency coun¬ 
terpart, the low-frequency viscosity is a non-perturbative 
quantity, and cannot be calculated by truncating the per¬ 
turbative series of Feynman diagrams to any finite order. 
It is thus necessary to sum an infinite set of diagrams, to 
all orders in a ee . This requires (i) that the Green’s func¬ 
tions are dressed by self-energy insertions and (ii) that 
vertex corrections are included. 

To bypass the tedious task of expanding the infinite se¬ 
ries of diagrams for the current-current response function 
to order q 2 , and then take the limits shown in Eq. (25), we 
use the definition of the shear viscosity given by Eq. (37), 
which we recall here: 


Vo = - Hm 

u ;—>-0 


3’rn[Xxy,xy(q = 0 ,^)] 

CO 


( 56 ) 
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FIG. 3. a) The diagrammatic representation of the current- 
current response function. The black filled circle represents 
the bare vertex A (0,£ d (we suppress the momentum-energy de¬ 
pendence for brevity), while the solid double lines are Green’s 
functions dressed by the self-energy. In the large-Alp limit it 
corresponds to the GW self-energy, which is depicted in panel 
b). Wavy lines represent the RPA screened interaction W. Fi¬ 
nally, the triangle represents the vertex function A 3 which is 
dressed by e-e interactions and satisfies the Bethe-Salpeter 
equation in panel c). Note that the form of the irreducible 
interaction I is uniquely determined by the choice of the self¬ 
energy, provided that A 3 satisfies the Ward identities 4 (see 
Fig. 4). 



FIG. 4. Feynman diagrams that contribute to the irreducible 
interaction I in Fig. 3. 


The greatest advantage of using Eq. (56) is that the 
hydrodynamic shear viscosity rjo can be extracted from 
the response function Xxy,xy{q,u) evaluated at q = 0 . 
Fig. 3a) shows the general expression of the stress-stress 
linear response function we calculate in what follows. 
The double solid lines represent Green’s functions dressed 
by the self-energy insertion shown in Fig. 3b). The wavy 
line represents the dynamically screened e-e interaction 
W(q,w). In the spirit of the large- Nf approximation, 
the self-energy can be approximated with its GW expres¬ 
sion. This, in turn, implies that the screened interaction 
in Fig. 3b) contains the usual RPA dynamical dielectric 
function e(q,uj). Herein lies our large- Nf approximation. 

The black filled circle on the left side of Fig. 3a) rep¬ 


resents the bare stress-tensor vertex 


= v f'~ 


k i S{ j Uk-,k + ) + k j SX>(k-,k+) 


(i) 


(57) 


To account for vertex corrections, we need to dress one of 
the two vertices in Fig. 3a). The dressed vertex is marked 
as “A” and it is dictated to satisfy the self-consistent 
Bethe-Salpeter equation represented in Fig. 3c). The 
choice of the quasiparticle self-energy [Fig. 3b)] and the 
requirement of fulfilling the Ward identities uniquely de¬ 
termine the irreducible interaction I. The diagrams that 
contribute to / are shown in Fig. 4. 

The present calculation closely follows that reported 
in Ref. 22, the main difference being that here we have 
stress-tensor vertices in lieu of current vertices. There¬ 
fore, the calculation of Ref. 22 should be adapted to the 
present case. These changes are explained in detail in 
Appendix E. Here, we briefly summarize the procedure 
we adopted to solve the problem posed by the diagrams 
in Figs. 3-4. 

As explained above, the choice of working with the 
stress-stress response function allows us to avoid the ex¬ 
pansion of the Bethe-Salpeter equation for the current 
vertex to order q 2 . We therefore set q = 0 from the 
very beginning of the present calculation. In what fol¬ 
lows we focus on the imaginary part of the stress-stress 
response function, which is the quantity that controls the 
low-frequency viscosity [recall Eq. (56)]. Taking the low- 
frequency and -temperature limits, it becomes evident 
that only the states in a narrow region of size ~ k^T 
around the Fermi surface give a significant contribution 
to the imaginary part of stress-stress response function. 
Moreover, in the limit ayr” 1 <C 2ep, whenever a prod¬ 
uct of two Green’s functions with the same momentum 
argument appears, one of the two must be considered as 
retarded and the other as advanced (schematically, we 
get products of the form G^G^). Terms containing 
products of two retarded or two advanced Green’s func¬ 
tions are responsible for subleading contributions in the 
limit cpr ee 1. Finally, in the spirit of the theory of 
normal Fermi liquids, each product G^G^ is approx¬ 
imated by a (5-function, i.e. the incoherent part of each 
Fermi-liquid Green’s function is neglected and we retain 
only its singular part. This corresponds to the so-called 
“quasiparticle approximation”, which is usually done in 
time-dependent perturbation theory to derive the Boltz¬ 
mann transport equation from the Keldysh equation for 
the Green’s function. Therefore, our theory describes the 
transport of momentum in the “semiclassical” regime, 
and neglects all quantum effects. 

The set of approximations described above allows us 
to dramatically simplify the expression of the Bethe- 
Salpeter equation. The latter can be then solved ana¬ 
lytically with the Ansatz? 7 (suppressing all the frequency 
and momentum dependencies) A = yA^ 0 ). In this way we 
require the dressed vertex to be proportional to the bare 
one. The Bethe-Salpeter equation becomes an algebraic 
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equation, which can be solved analytically. We get 

u + i/T ee 

7" = —rr7— • 

U) + l/T v 

The expression for r v reads (see also Appendix E) 


(58) 


1 

T v 


/»+ oo 


W 

Vn 


dui 


1 - n F (^fc,+ - op 
1 - exp(-/3w) 7 0 


(*+oo 


dq q 


e(q,u>,T) 


3m[x (0) (q, W, T)] A ++ (/c, q, u) 


x 4 

where 

A ++ = 


1 - 


q 2 - cc 2 /ug 
4k(k — tu/v p) 


q 2 ~ u 2 /v 2 
2 k(k — oj/v-p) 


(59) 


4(k — oj/vp) 


v Fy /\(2k - ui/vp) 2 - q 2 ](q 2 - w 2 /uf) 


1 - 


9 2 - U] 2 /l 


4 k(k — co/vp) 
x ©| [(2k - w/u F ) 2 - q 2 ] (q 2 - w 2 /u|)} . 


(60) 


Eq. (59) describes the scattering of a quasiparticle with 
momentum k = kpx and energy equal to the Fermi en¬ 
ergy into a quasiparticle with momentum \k + q\ = kp. 
In this process the whole Fermi liquid is excited. The 
many-particle excitations created during the scattering 
event are encoded in the imaginary part of the density- 
density response function 3m[xin(g,w,T)]. The last 
line of Eq. (59) can be shown to be proportional to 
1 — cos 2 ((fik +q ) [see Appendix E —to get Eq. (59) an inte¬ 
gration over the angle of q has been performed]. At low 
temperature and in the limit of a ee —> 0 Eq. (59) gives 

1 „ 87t (kpT) 2 

3 hsp 


' a ee ln(a ee ) • 


(61) 


Note that the matrix element 1 — cos 2 (ipk+ q ) vanishes 
when k + q is either parallel or antiparallel to k (recall 
that k is along the x direction), and it is maximum for 
transverse excitations, for which k + q is perpendicular to 
k. Therefore it suppresses the contribution coming from 
the region of small transferred momenta q ~ 0, which 
is responsible for the logarithmic-in-temperatures correc¬ 
tion to the quasiparticle lifetime. Indeed, at low temper¬ 
ature 1 /T ee oc — T 2 ln(T), with a coefficient of proportion¬ 
ality which is independent of the coupling constant a e e 
(see Ref. 5 for more details). Therefore, because of the 
presence of the matrix element 1— cos 2 (<pk+ q ), l/r v oc T 2 , 
and the coefficient of proportionality depends on a ee . 

Note also that, since the derivation of the viscosity 
transport time does not rely on the linear band dispersion 
of MDFs, a similar conclusion can be drawn for a Galilean 
invariant parabolic-band 2D electron gas. Also in the lat¬ 
ter case there is no logarithmic-in-temperature correction 
to the low-temperature behavior of r v . Indeed, the ma¬ 
trix element suppresses both the regions q ~ 0 (forward 
scattering) and q ~ 2 kp (perfect backscattering), which 
are responsible for the aforementioned correction 4 . 



FIG. 5. (Color online) Panel a) shows the low-frequency kine¬ 
matic viscosity no = Vo/(nm c ) of a 2D MDF liquid (in units 
of m 2 /s) as a function of the excess carrier density n (in units 
of 10 12 cm -2 ). The three curves refer to different values of 
the temperature, i.e. T = 150, 200, and 270 K. In this plot 
cc e e = 0.5. Note that in the range of densities shown in this 
plot the temperature is always smaller than Tp = ep/kp. 
This, in turn, implies that the low-temperature approxima¬ 
tion (T -C Tp) for 1/t v in Eq. (59) is valid. Indeed, the 
minimum Fermi temperature in this plot, corresponding to 
n = 0.2 x 10 12 cm" -2 , is Tp ~ 605 K. Therefore, in this plot 
T/Tp < 0.45. Panel b) The low-frequency kinematic viscosity 
is shown as a function of temperature T (in units of K) for 
an excess carrier density n = 0.4 x 10 12 cm -2 (corresponding 
to Tp ~ 855 K). Inset: a logarithmic plot of no in the same 
range of temperatures as in the main panel. Note that no 
grows like T~ 2 for T <C Tp. A crossover to a different power 
law, ~ T ~ 3 / 2 , is evident as temperature increases beyond the 
T <C Tf regime. 


V. RESULTS AND CONCLUSIONS 


In Fig. 5a) we illustrate the carrier density dependence 
of the low-frequency kinematic viscosity n 0 = rj 0 /(nm c ) 
of a 2D MDF liquid in a doped graphene sheet. We 
show three curves corresponding to different values of 
temperature T, i.e. T = 150, 200, and 270 K. Data 
in this plot have been obtained by setting a ee = 0.5. 
In Fig. 5b) we plot the same quantity but viewed as a 
function of temperature T (in units of K), for an excess 
carrier density n = 0.4 x 10 12 cm -2 . In the inset we also 
show a logarithmic plot of n 0 , and we compare it with 
the power laws T~ 2 and T -3 / 2 . While the kinematic 
viscosity goes as T~ 2 for very small temperatures, we 
find that it goes as T -3 / 2 for the larger temperatures 
shown in Fig. 5b). Note that the kinematic viscosity of 
the 2D MDF liquid in graphene is much higher than that 
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FIG. 6. (Color online) Panel a) shows the kinematic viscosity 
i/ w = tj u /(nm c ) (in units of m 2 /s) of the 2D MDF liquid 
in a doped graphene sheet—Eq. (1)—as a function of w (in 
units of the Fermi energy £f). In this plot we set n = 0.2 x 
10 12 cm -2 , a ee = 0.5, and we show results for three different 
temperatures: T = 150, 200, and 270 K. All the curves tend 
to the finite (albeit small) value at large t o, whose magnitude 
can be extracted from Fig. 2. In this plot we subtracted 
the “conductivity term” proportional to BN c (ct eB ) from the 
definition of r/oa given in Eq. (53). Panel b) Same as in panel 
a) but for n = 0.5 x 10 12 cm -2 . Panel c) Same as in panels 
a) and b) but for n = 1.0 x 10 12 cm -2 . 


of classical fluids 1 . 

In Fig. 6 we illustrate the frequency dependence of the 
kinematic viscosity = r]^/ ( nm c ) of the 2D MDF liquid 
in a doped graphene sheet. The quantity rj w is defined 
in Eq. (1). We remind the reader that w here represents 
the frequency of the external perturbation. In Fig. 6 a) 
we fix the excess carrier density n = 0.2 x 10 12 cm -2 , the 
coupling constant a ee = 0.5, and we show three curves 
for T = 150, 200, and 270 K. In Fig. 6 b) and c) we fix 
instead the carrier density at n = 0.5 x 10 12 cm ” 2 and 
n = 10 12 cm -2 . Fig. 6 shows that viscosity corrections to 
the lifetime of plasmons in the high-frequency (e.g. mid 
infrared) regime are totally negligible. However, when 
the plasmon frequency is in the Terahertz (THz) spec¬ 


tral region [cj/ef <0.1 in Figs. 6 a)-c)] the viscosity of 
the electron liquid leads to corrections 38 to the plasmon 
lifetime that may be comparable to those due to electron- 
impurity 19 and electron-phonon scattering 20,39 . There¬ 
fore, a careful comparison between measurements of the 
lifetime of THz plasmons in high-quality graphene sam¬ 
ples and theoretical predictions can be used to extract 
the value of v^. 

In summary, we have calculated (i) the high-frequency 
bulk and shear viscosities—Eq. (53)— and (ii) the vis¬ 
cosity transport time—Eq. (59)— of the two-dimensional 
massless Dirac fermion liquid in a doped graphene sheet, 
as solely due to electron-electron interactions. As ex¬ 
pected, the bulk viscosity vanishes. The shear viscosity 
is instead finite and is proportional to the excess carrier 
density of the doped system. Note that, since the bulk 
viscosity vanishes, the high-frequency shear viscosity 77 ^ 
can be directly estimated by measuring the lifetime of 
Dirac plasmons 14 in high-quality encapsulated graphene 
sheets 39 , by carrying out similar experiments in the Ter¬ 
ahertz spectral range. Finally, from the knowledge of the 
high-frequency values of the shear viscosity and modu¬ 
lus, and of the viscosity transport time, we extracted the 
shear viscosity r at all frequencies, using the interpola¬ 
tion formula given by Eq. (1) and derived in Appendix B. 

The low-frequency shear viscosity is equal to the high- 
frequency shear modulus Soo = ne F /4 multiplied by a 
“viscosity transport time” r v , which we microscopically 
calculated in this Article. We showed that r v is both 
quantitatively and qualitatively different from the quasi¬ 
particle lifetime 5,6 r ee . Indeed, in the low temperature 
limit 1/ree oc —T 2 ln(T/Tp), and the coefficient of pro¬ 
portionality is independent of the coupling constant a ee 
(as shown in detail in Ref. 5). Conversely, we proved that 
1/t v ocT 2 and that it depends on a ee . Transverse excita¬ 
tions, in which a quasiparticle is scattered at 90° with re¬ 
spect to its initial direction of motion, are responsible for 
the dominant contribution to the (shear) viscosity trans¬ 
port time. This is reflected by a matrix element in the 
expression of l/r v , which suppresses the contributions 
due to both the forward scattering and perfect backscat- 
tering of the quasiparticle [i.e. the perfectly longitudi¬ 
nal channels). The same matrix element is expected to 
show up in the expression for the viscosity transport time 
of a Galilean invariant (parabolic band) two-dimensional 
electron gas. Therefore, also in the latter case we expect 
1 /t v oc T 2 for T < T f . 


VI. ACKNOWLEDGEMENTS 

A.P. and G.V. were supported by the BES Grant DE- 
FG02-05ER46203. M.C. acknowledges support by MIUR 
through the program “FIRB - Futuro in Ricerca 2012” 
- Project HybridNanoDev (Grant No. RBFR1236VV). 
M.P. is supported by MIUR through the program “Pro- 
getti Premiali 2012” - Project ABNANOTECH and by 
the European Community under the Graphene Flagship 








14 


program (contract no. CNECT-ICT-604391). 


Appendix A: The equivalence of hydrodynamic and 
microscopic approaches in the Fermi liquid regime 


In this Appendix, we demonstrate the equivalence of 
hydrodynamic and microscopic approaches for the 2D 
MDF liquid in a doped graphene sheet, in the low- 
temperature Fermi liquid regime. 

To this aim, let us focus on the non-interacting equa¬ 
tion of motion for the current operator 

i&tjq.i — 

= -Vp qiflq 

— 2ivp E ^Lq/2,a( fc X °'a/3)V’fc+q/2,/3 • 

k,<x(3 

(Al) 


We remind the reader that we are interested in calcu¬ 
lating the current-current response functions to order q 2 , 
in order to take the limits of Eq. (25). This in turn 
implies that we can retain only the term of order q of 
Eq. (Al), and replace the first term on its right-hand 
side with qifi q -A qifi q= o- Since n q= o is the total num¬ 
ber of particle, it is exactly conserved by any microscopic 
process, and therefore does not contribute to any linear 
response function. Thus, we neglect it in what follows. 
We now use the fact that we are interested only in the 
states in conduction band around the Fermi surface. Af¬ 
ter a rotation of the second term on the right-hand side 
of Eq. (Al) to the chiral basis, we retain only the states 
that have chiral indices A = A' = +. We thus get 


dtjq.i = - 2 4E4-g/2,+ ( fc X ^)Cfc+q/2,+ 

k 

x5<t(fc-<l/2,Hg/2), (A2) 

where S^,(k,k') is defined in Eq. (13). Therefore, 


<S++(fe — q/ 2, k + q/2) = isin 


@k—q/2 @k-\-q/2 


1 q x k 

2 k 2 


+ 0(q 2 ). (A3) 


Here we used that, close to the RT-point of the Brillouin 
zone 

k v 


Vfcdfc = Vj 


arctan | - 


ky k x 

W'W 


(A4) 


Plugging Eq. (A3) into Eq. (A2), after some simple alge¬ 
bra we can rewrite it to 0 (q 2 ) as 


'i^t jq.l — E] Ofe.H 




v F qj 


k 2 Jc 1c 

"'x rx, x rx "y 


m c k F \ k x k v k 2 


x^y ryy 


O'* 


v F qi 


EC 


j/2,+ Cfc +q/ 2 >+ 


(A5) 


The term on the second line is proportional to the to¬ 
tal number of particles in the conduction band. At low 
frequencies, since interband processes are strongly sup¬ 
pressed, the number of particles in each band is con¬ 
served. Therefore, such term does not contribute to the 
response functions and can be neglected. On the other 
hand, we note that the matrix elements of the current, 
when projected at the Fermi surface, read «S++(fe, ft) = 
ki/k f- Therefore Eq. (A5) can be rewritten as 


i&tjq,i 


E 


C fe , + Cfe,+ 


v F gj kiS%l(k, ft) + kjSP+jk, k) 
m c 2 


(A6) 


Since (i) we always consider states around the Fermi sur¬ 
face, (ii) we do not allow interband transitions, and (iii) 
we use the equation of motion to order q, after a rotation 
back to the pseudospin basis we can rewrite Eq. (A6) as 


i^tJlq. / 


m c 


k,a,/3 




q/2,u 


k iKp + k 0<P 


4’k+q/2,0 ■ 

(A7) 


The quantity on the right-hand side is exactly the bare 
stress tensor operator of Eq. (31). Using Eq. (A7) we can 
derive the analogous of Eq. (39) for an MDF liquid and 
prove the equivalence of the two approaches [although in 
the approximate sense explained after Eq. (A6)]. 


Appendix B: The generalized relaxation time 
approximation 

The Generalized Relaxation Time Approximation 
(GRTA) offers a simple and powerful way to interpolate 
between the collisionless and the collisional regimes of 
linear response functions. This theory has its roots in 
Mermin’s work on the density-density response function 
of the electron gas 40 , but it is applicable to a broader class 
of response functions and problems. The classic deriva¬ 
tion is from the solution of the Boltzmann equation with 
a collision integral that describes relaxation towards a 
local equilibrium distribution function. 

A local equilibrium distribution function (LEDF) has 
the same form as the equilibrium distribution function (a 
function of constants of the motion) but is evaluated at 
the local values of the densities of the conserved quan¬ 
tities. Thus, in the presence of impurity scattering, the 
LEDF depends on a local chemical potential and a lo¬ 
cal temperature, determined by the local particle density 
and energy density respectively. Momentum is not con¬ 
served, and therefore the LEDF has zero average momen¬ 
tum. This is the situation considered in Mermin’s orig¬ 
inal paper. On the other hand, if impurity scattering is 
absent or negligible on the time scale of particle-particle 
collisions, then the LEDF depends also on a local drift 
velocity, determined by the local momentum density. 
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Regardless of whether electron-impurity collisions or 
electron-electron interactions dominate (we assume that 
it is always either one or the other) one can distinguish 
between (i) a collisionless regime in which the frequency 
of the macroscopic motion is much higher than the col¬ 
lision rate, and (ii) a collisional regime in which the fre¬ 
quency of the macroscopic motion is much lower than the 
collision rate. 

When electron-impurity collisions dominate, the colli¬ 
sional regime is known as diffusive regime: there is a di¬ 
rect coupling between the current and the gradient of the 
density. Similarly, when electron-electron collisions dom¬ 
inate, the collisional regime is known as hydrodynamic 
regime: there is a direct connection between the stress 
tensor and the gradient of the velocity field. 

In the diffusive regime a drift velocity can only arise 
from the deviation of the distribution function from the 
LEDF, whereas in the hydrodynamic regime a drift veloc¬ 
ity is already present in the LEDF. The diffusive regime is 
intrinsically limited to drift velocities that are small rela¬ 
tive to the Fermi velocity - it is basically a linear theory 
of drift. Whereas the hydrodynamic regime is suitable to 
describe large drift velocities, which cannot be treated in 
the linear approximation. 

Let us formulate the general theory of the GRTA for 
a generic response function x(<Z, w). We define the relax¬ 
ation function K(q,ui) as follows 41 : 


interactions only to the extent of renormalizing the Fermi 
liquid parameters. Thus the quasiparticle lifetime and 
all the transport times, whether they be due to electron- 
electron, electron impurity, or electron-phonon collisions, 
are assumed to be infinite: we are in the collisionless 
regime. The associated “high-frequency” relaxation func¬ 
tion is denoted by K ao (q,uj). We posit that the relation 
between the full interacting relaxation function and the 
“high-frequency” one has the following form: 


1 

K(q,u) 


1 

Kqo (q,u + y) 


V T {q,w) 


(B4) 


where the replacement ui —> to + i/r accounts for self¬ 
energy insertions, and V T (q,oj) for vertex corrections. 
The two corrections 1/r and V T are not independent of 
each other, however. They must be chosen consistently, 
in order to satisfy the conservation laws. In the following 
subsections we consider the case of the density response 
in both disorder and clean systems. Therefore, all the 
response functions that appear in what follows should 
be regarded as density-density ones (we suppress their 
indices for brevity). 


1. GRTA in a disordered system 


X( 9 ,w) = x(g,0) [1 +iuK(q,uj)} (Bl) 


Since we will always be working at small q (q <^. ftp), the 
g-dependence of %(g, 0) will be ignored: x(<?,0) ~ %o- 
Here \ 0 is (minus) the density of state of quasiparti¬ 
cles, and is thus renormalized by electron-electron inter¬ 
actions. The relaxation function can further be expressed 
in terms of a generalized frequency- and wave vector- 
dependent relaxation time T(q,u>) for the quantity under 
consideration. Its general structure is 

K(q,u)= . * i (B2) 

IU) + 


This form guarantees that the finite frequency response 
function y(g,w) vanishes if the generalized relaxation 
time T is infinite, i.e., if the quantity under consideration 
is conserved. From Eqs. (B1)-(B2) it is indeed immediate 
to get 


xOlw) = 


Xo 


1 — ivjT(q, ui) 


(B3) 


Let us first apply the theory developed so far to a sys¬ 
tem in which the dominant mechanism is the electron- 
impurity scattering. Since the only conserved quantity is 
the particle number, we now consider the relaxation func¬ 
tion of density fluctuations. At q = 0 we have that the 
non-interacting relaxation function is (0, ui) = — iu.i. 
Then we simply choose V T (q,oj) = 1/t so that Eq. (B4) 
yields K~ l { 0,w) = — iu, as it should be due to particle- 
number conservation. The inverse density-density re¬ 
sponse function, as calculated from Eq. (Bl), reads 

1 = 1 iu K(q,u) 

x(?>w) Xo Xo 1 + iwK(q,u) 

Substituting in Eq. (B5) the Ansatz (B4), we get 

1 _ 1 iu} K 00 (q,u} + f) 

X(q,u) xo Xo l + i[u + iV T (q,u)\(q,u + i) 

(B6) 


It is essential to realize that T(q,ui) is the lifetime of a 
collective mode (say a density fluctuation) and as such 
is conceptually different from the microscopic scattering 
times (quasiparticle lifetime and transport lifetime) upon 
which it may depend. In fact, we will see that the form of 
T{q,(jS) depends crucially on the nature of the response 
function under study. 

Let us now introduce the “high-frequency” response 
function \oo (g,w). The latter includes electron-electron 


At the same time, taking the limit r —► oo in Eq. (B5), 
and replacing ui —> u> + i/r we have 

1 _ 1 
Xoo (<7,w+y) Xo 

*(w+y) Koo (g, UJ + y) 

Xo 1 + * (w + y) Aoo (g, U; + y) 

(B7) 
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Combining Eqs. (B6) and (B7) we finally get 


2. GRTA in a clean system 


L +_H_ ( _x°__ _ 1^1 

io (w + y) Xo \ Xoo ( q , w + y) ) 

r{q,u) - y] / Xo _ A 

( w + t) \X°° (^ + t) ) 

(B8) 

Setting now V T (q,uj) = 1/r we get 

1 _ £ 1 w 1 

X'HH (w + y) Xo (w + y) Xoo (q, w + y) 

(B9) 


1 

xfe w) 



We now discuss the case of a clean Galilean-invariant 
system in the so-called “hydrodynamic regime”, in which 
the main mechanism of quasiparticle relaxation is given 
by electron-electron interactions. Replacing r —> r ee , 
from Eqs. (B9) we get the “naive” diffusion constant 


D T = 


Pw 1 

* oi (" + i) 


(B14) 


However, since electron-electron collisions conserve the 
total momentum the true diffusion constant must be di¬ 
vergent in the limit of ui —> 0, i.e. it must be 


which is precisely the Mermin’s formula for the density- 
density response function 40 . In a Galilean-invariant elec¬ 
tron gas in the diffusive regime 


D t r = 


Xo ’ 


(B15) 


1 

x(q,u) 



(BIO) 


where D T = —'D w t/(x o) is the diffusion constant, V w 
is the renormalized Drude weight ( n/m for the two- 
dimensional electron gas or kpvp/fi for massless Dirac 
fermions), and we have used the oj —>■ 0 diffusive form of 
the density-density response function Xoo (?)W + i/r) ~ 
-V w q 2 r 2 . 

However, it is well known that the relaxation time 
appearing in the diffusion constant, commonly denoted 
by r tr , is different from the self-energy relaxation time 
r (quasiparticle lifetime). The conservation of particle 
number, as implemented in the q = 0 vertex correction, 
is not sufficient to produce the correct transport relax¬ 
ation time. To capture the difference between r and T tr 
we must allow for the ^-dependence of the vertex correc¬ 
tion in our Ansatz (B4). Namely, we set 

V T (q,uj) = - + D T q 2 . (Bll) 

T 


The coefficient D T has the dimensions of a diffusion con¬ 
stant and its value is determined by the requirement that 
the final formula for x(9, uj) has a diffusive limit with the 
correct relaxation time, i.e., r tr . With the choice (Bll) 
for the vertex correction our Eq. (B8) gives 


1 1 iui 1 

XHH “ Xo Xo ( D t - D T )q 2 


(B12) 


It is evident from this result that we obtain the correct 
diffusive limit if and only if 

D T =D T (l- (B13) 


Thus, the coefficient D T in the vertex correction is simply 
the difference between the “naive” diffusion constant and 
the “true” diffusion constant, in which r is replaced by 
the transport lifetime r tr . 


Therefore the role of D r is played by 

D r = D r -Ar = -— - 7 1 . v • (B16) 

X'O r ee w ( u) + 

Substituting this in Eq. (Bll) for the g-dependent ver¬ 
tex correction, and then putting this vertex correction in 
Eq. (B8) we finally get 


ui 


» ("+£) W («■"+£) 

iV w q 2 / (wT ee ) 


J- Xo 


X < 1 + 


Hi)' 


(<^ +a) 


1 

Xo 


(B17) 


We now introduce 1/Xoo (<?,w) to be the q° coefficient of 
the expansion of 1/Xoo (<7, H> ie. 


w 


Xoo (g,w) A<? 2 Xoo (q,u) 

This allows us to rewrite Eq. (B17) as 


(B18) 


1 


1 


ui 


x(q,v) xo + ^ 
l 


l - 


Xo 


v Xoo(</,w + £) *»(9 > w + £). 


X U- 


U) 


1 + 


WTeeXO 


Xoo + 


(B19) 
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Expanding Eq. (B19) for q <C ui, l/r ee £f, and retain¬ 
ing all the terms up to the order q°, we get 


1 

Xfew) 




i 

WT ee X0 

(B20) 


To make a clearer connection with the solution of the 
Boltzmann equation in the relaxation time approxima¬ 
tion we rewrite 



OJ ioj P w g 2 

oj + ^ T ee V w q 2 (w + */r ee ) 2 

(B21) 

The first term reproduces the Mermin formula, and the 
second term is the correction required to ensure momen¬ 
tum conservation, i.e. 



(B22) 


To extract the hydrodynamic coefficients from 
Eq. (B22) we now observe that 


1 w 2 

b~ + ( 2 2 ) 

Soo (cj) 

Xoo (q,v) 2? w q 2 

n 2 \ d) 

n 2 


(B23) 

where the bulk modulus Boo is real, frequency- 
independent, and related to xo by the well-known com¬ 
pressibility sum rule 


Xo = 


n 

K 


(B24) 


while Soo(u}) = Soo — iurioo is the complex shear modulus. 
Plugging Eqs. (B23) and (B24) into Eq. (B22) we get, 
after some simple algebra, 

1 _ cu 2 

X'(<?,w) “ V w q 2 

Boo /_ ^ \ ^ \ ^oo 4“ i/^~ee) 

n 2 \ d) y w + J n 2 

(B25) 

These result does not matches with what we found in the 
main text of the paper. Indeed, the correct result should 
be 


X(q, w) V w q 2 

Bqq _ 2 \ u Sqq(u}) 

n 2 \ d J oj + i/r v n 2 

(B26) 


To get the correct result it is necessary to take into ac¬ 
count the corrections to order g 4 to the vertex correction 
V T (q,oj). We therefore set 


V T (q,uj) 


1 

T 


D r q 2 


i^(q,uj) D 2 g 4 
Xo (w + i/r ee ) 3 


,(B27) 


where S(g,w) is to be determined. Plugging Eq. (B27) 
into Eq. (B8) and expanding it up to order g , after some 
straightforward algebra we get 


1 


OJ 


X(<7,w) 

Bo 




2 - 



Boo (^0 4“ ^/^ee) 



(B28) 


It is clear that it is possible to choose ^(q,oj) in such a 
way as to match Eq. (B28) with Eq. (B26). 


Appendix C: Calculation of F\ and j\, q 


The operator Fi was calculated in Ref. 18 and reads 

1 
2 


A,A ' 


X C fe_,A Cfc +A' C fc' 


(Cl) 


where 


,q’) = 


X>aa '(fc-,fc+)V(^,fc-) 

Ek-,\ £fe+,A' H” £k' + ,fi £k'_,n' 

~ (C2) 
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To simplify the notation we have introduced k± = k ± 
q '/2 and k' ± =k'± q 1 /2. 

The transverse component of the operator j-\ q is ob¬ 
tained from its definition of Eq. (41). After some manip¬ 
ulations, along the lines of Ref. 18, we get 


jl,q,T 



-y-( T ) " 


VT,, 


(T) 


(C3) 


where 


'Y'(T) _ 

qq' 


E a , q)4_ - q/ 2 ,\Ck + +q/ 2 ,X' 
fc, A,A' 


(C4) 


with 


M x ,\'(k,q',q) = 

E 


w + £ fc_+q/2,p — £ k^~q/ 2,A 

(C5) 


2? 


Ap 


k --2’ k +-2) s %( k+ -r k + 


S$ ( 


w + e fc++q/2,A' — £ k + ~q/2,p 

k 


k --% k ~ + \) V P* ( 


9 , 9 

2’ fc++ 2 


We stress that Eqs. (C3)-(C5) are valid only for their use 
in Eq. (42), and become exact in the low-energy MDF 
continuum limit, which will be taken momentarily. The 
same limit restores rotational invariance: we therefore 
have the freedom of fixing the direction of the wave vector 
q arbitrarily. In deriving these equations we have taken, 
without any loss of generality, q = qy which implies that 

Jl,q,T = jl.q,x- 

Following the same steps outlined in Ref. 18 it is possi¬ 
ble to show that the operator in Eq. (C4) can be written 


as (the details of this derivation are given in App. D) 

t (T) =V f v F g \ q x q y n f r -l' 2 \ 

qq ' a h,y\“ 2k ? I 9' 2 V 4 k\) 

X {q x Sa : y + 9y<5a,y)] + ^,3 J jq',a 

= Y, Fi T) (9,9')jq',a ■ (C6) 

a=x,y 


The main differences between Eq. (C6) and the corre¬ 
sponding expression that can be found for a 2DEG 36 are: 
(i) the factor 1 — q' 2 /(Akp), which is due to the chirality 
of the MDF eigenstates and suppresses backscattering at 
the Fermi surface 18 , and (ii) the last term in curly brack¬ 
ets, which is finite even in the long-wavelength q —> 0 
limit 18 . 

For a sake of completeness we recall 18 that the first- 
order correction to the longitudinal current operator 
is formally identical to Eq. (C3), with the longitudi¬ 
nal counterpart [T^,] of the operator T^, taking 

its place. We recall that T*l L L is also proportional 
to the untransformed current operator, i.e. ± q q ' = 

Ea= Il9 r « L) (9.9')j,', a , where 

/v_^F9x 

’ q> W 2 [q' 2 k F 
q ' 2 A 

+ 4 v F k* a ’ x 

After the change of variables q' —► —q' in the second 
term on the right-hand side of Eq. (C3), the latter can be 

rewritten as q ■ j\ q = J2 q ' v q'^ qq ' ,£l -q'■ A major sim- 
plihcation is suggested by the analysis of the Feynman 
graphs contributing to the non-interacting spectrum of 
9 ■ ji.q- As shown in Ref. 18 the disconnected graphs 
contain two independent sums over the number Nf of 
fermion flavors, whereas the connected ones contain only 
one such sum. We conclude that the disconnected graphs 
dominate in the large-W limit. The final formula for 
the two components (longitudinal and transverse) of the 
current-current response function, which is exact to sec¬ 
ond order in e-e interactions and in the large -Nf limit, 
is 




Sm[xf(g,w)] = - 


a.,(3=x,y 1 


d 2 q' 

(2tt) 


2< 


' J o ^ { r a } (<?’ q') 1 ^ (-9, -q')3rn[xW(q',u')]5tm[x { y jf) (q',uj- a/)] 


+ r V\q,q')r ( g\-q,q')%m[x { °] C '{-q',u>')} (q', u - a/)]} , 


(C8) 


In this equation x™(g,w), (q,w), and Xnf (9i w ) and density-current response functions of a 2D gas of 

JOlJp * J OL , p\ 

are the non-interacting density-density, current-current, MDFs. The quantities {r„ (q, q');a = x,y;£ = L, T} 
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are defined in Eq. (C6) and (C7). We stress that the 
imaginary parts of the three linear-response functions 
Xnl(q,u), Xj° ) i/3 (<Lw), and XnjS 1 ?,<*>) do not depend 
on any ultraviolet cut-off in the low-energy MDF limit. 
Moreover, since in the limit of ta — > 0 the integral over 
q' is naturally restricted to 0 < q' < 2fc F , no ultraviolet 
regularization is needed in Eq. (C8). The only pathology 
of the integral in Eq. (C8) appears in the infrared q' —0 
limit, due to the 1/q' singularity of the Coulomb poten¬ 
tial v q '. This problem is cured by screening, as discussed 
in the main text. 


Appendix D: Details of the manipulation of T^, 

In this appendix we approximate the expression of 
the operator in Eqs. (C4)-(C5) by taking the limit 
vpq <C w <C 2e F . We will try to slowly guide the reader 
through the many steps of this lengthy process. 

To begin with, in the long-wavelength q —0 limit we 


can write 


1 

w + £ k±+q/ 2,A — e k±-q/2,p 


"A ,p 


1 

W 


+ (1 - ^A.p) 


q de k± 
u> 2 dk y 

1 

La -f- 2 As‘f 


+ 0(q 2 ) . (Dl) 


We then observe that in the regime of interest in this 
Paper, i.e. v F q <C oa <C 2e F , the particle-hole states 
created by the operator T^, are energetically close to 
the Fermi energy. The band indices A, A' on the right- 
hand side of Eq. (C5) are therefore constrained to take 
the values A = A' = +1 (recall that £ F > 0). 

Note that the “virtual state” p , over which the sum 
on the right-hand side of Eq. (C5) runs, can be either 
in conduction (p = +1) or valence (p = —1) band, even 
though the states labeled by the band indices A and A' 
are bound to the Fermi surface. 

We first simplify Eqs. (C4)-(C5) by using Eq. (Dl). 
We are naturally led to define 


M intra (k,q',q) = M ++ (k,q',q) | p=+1 = cos 


?fc-q/2 — 0k + ~q/2 


— COS 


fc+qr/2 t 'fc + +qr/2 


cos(0 fc J 


1 v • (a \ 

-y Sin(0fc_ ) 

r ,1 t, 


cos(6»fc ) - - ^ sin(0 fc ) 

\_UJ W“ 

0(q 2 ) , 


and 


M inter (fe, q\ q) = M ++ (k, q 1 , q)\ p= _ 1 = sin ^ 


2e F 


■ Sill 


^fc+q/2 @k + + q/ 2 


Sill 


0k-—q/2 ^k+—qj2 


0 k --q/2 + 0 k _+q/2 


Sill 


®k + -q/2 + 0k + +q/2 


+ 0(q 2 ) 


(D2) 


(D3) 


so that in the limit v F q -C to -C 2£ F we have 

M ++ (k,q',q) = M intra (fc, q', q) +M inter (fe, q', q) . (D4) 

In writing Eq. (D3) we have taken the limit ua —> 0 in the 
second term on the right-hand side of Eq. (Dl). More¬ 
over, in obtaining Eq. (D2) we have used that 


cos 


+ = cos{ekJ + 0[tf) (D5) 


and 


de , 


k± ,A 


dk„ 


Au F sin (0 k± ) . 


(D6) 


The last equation becomes exact for k close to the K 
point of the BZ and therefore in the low-energy MDF 
limit. 

Clearly we can carry out further approximations, rely¬ 
ing on the fact that we are interested in the low-energy 
MDF limit. Eq. (D2) can be further simplified by noting 


that 


cos 


0 k _±q/2 — 0k + ±q/2 \ 
- 2 - 1=1 


COS 


'k——q /2 ^fc|+qr/2 


Q . 

--sm 


Ok^ - Ok, \ dOk 


dk. 


+0(q 2 ) , 

which leads to 

COs(dfc^) COS [0 k _ ) (0k_ — q/2 0kj r -\-q/2 


V 

(D7) 


-Mintra. — 


v F q 


■ COS 


aa 


—y [cos(6» fc _) sin(6» fc _) - cos (0 k+ ) sin(6» fc+ )] 


LJ 


X COS 


Ok- - Ok + 


q <9[sin(6» fc _) - sm(0 fc J] . (0 k _ - d k 

-sm 1 


2c o 

0(q 2 ). 


dk T 


(D8) 






























In the first term on the right-hand side of Eq. (D8) we 
can approximate 
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cos(0 fe+ ) - cos(0 fc _) 


kp 


(D9) 


while the second term on the right-hand side of Eq. (D8) 
becomes 


[cos(0 fc _) sin(0 fc _) - cos { 0 k+ ) sin(0 fc+ )] cos 


V k — 8 k 


k x Q y + k y q x 


4 


cos 


0 k ~ Ok. 


q'y COS (0 k ) + cos (0 k _) q' sin(0 fc+ ) + sin(0 fc _)' 


xcos 


fcp 2 

- Ok 
2 

„/2 


=~ 1 - 




fcp 


sin 


Ok + 0fc+ 


&F 


COS 


0fc + 0fc 4 


(DIO) 


Finally, the derivative in the third term on the right-hand 
side of Eq. (D8) reduces to 


Again, we used the fact that j q = VF& q close to the K 
point of the BZ. 

Summing Eqs. (D12) and (D14) we finally get 


-( T ) _ v F q x 

1 r, — t ll q+q' 


q,q 


fcp a; 


'Y'/ 


qq' 


(D15) 


with 


Y'/ _ 


vf q 
w 2 
„/2 


4uf fcp 


1 - 


Jq' ,x 


4fcp 


jq'.y + u Jq'.x 


(D16) 


The first term on the right-hand side of Eq. (D15) can 
be further manipulated. Indeed, when it is introduced in 
Eq. (C3) it gives 


2wfcp 


V 


'x n q+q' n —q' qx n q' n q-q' 


Vf ST' r / / i 

7^ 2 n q+q' n -q' V ^ 
q' 


vf q 


q.j; q„ 


(D17) 

q' 


d[sin(6> fc _) - sin(0 fc _)] 
dk y 


d(q'y/k F ) 
dky 


(Dll) 


Introducing Eq. (D8), approximated according to 
Eqs. (D9)-(D11), back into Eq. (C4) we get the “intra- 
band” contribution to the operator which reads 


T 


(T,intra)_ 

q,q' 


VFq x 


vf q 


k F co q+q ' cc 2 



X 


Jq'.y 



+ 0(q 2 ) . (D12) 


We remind the reader that, without loss of generality, we 
have taken q = qy. Here we used that j q = uf<t q close 
to the K point of the BZ. 

We now consider Eq. (D3). Steps similar to what sum¬ 
marized above for the intra-band contribution to T^, 
yield 


-Winter = x—[sin(0 fc _) - sin(0 fc )] sin 
Z£f 


o k - o k . 


i 

= — sm 
£f 
„/2 


Ok — Ok, 


cos 


Ok 


2 

0 k+ 


q ~ S$(fc_,fc+). 


4zipfc| 


(D13) 


Once Eq. (D13) is introduced into Eq. (C4), it gives the 
“inter-band” contribution to the operator Y^,, i.e. 


Here, we performed the shift q' —► q + q' in the term 
proportional to h q >h q - q ', using that q = qy, and we took 
the small-q limit in the last line of Eq. (D17). Finally, 
using the continuity equation 

ujfi q 'Ti— q ' — q • j q 'Ti— q ' -I - n q >q • J— q ' . (D18) 
we can rewrite Eq. (C4) in the form of Eq. (C6). 


Appendix E: The viscosity transport time 


We now guide the reader through the all-order di¬ 
agrammatic calculation of the low-frequency viscosity, 
from which we extract the corresponding transport time 
r v . We remind the reader that the low-frequency shear 
viscosity is given by 


Vo = - lim 

uj —>-0 


%rn[x X y, X y(q = 0,<u)] 

UJ 


(El) 


where w) is the stress-stress response function. 

Figs. 3-4 summarize the all-order microscopic calculation 
of Xot 0 ,^(q.u), which is given by [Fig. 3a)] 


k,e n ,A,A' 

x A x’y 3 ) {k~,k + )Gy(k + ,i£„+iuJm) 
x A ^\k+,i£ n + iu mi k— i i£ n ) • (E2) 


-^(T, inter) 

1 q,q' 


4^^ + ° ( ’ 2) - 


(D14) 


Here G\(k,i £) is the Green’s function (represented by a 
double solid line in Figs. 3-4) dressed by the self-energy 
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insertion of Fig. 3b), k± = k±q/ 2, e n (w m ) is a fermionic 
(bosonic) Matsubara frequency, and A and \' are band 
indices. In Eq. (E2) we defined the bare vertex (repre¬ 
sented as a solid dot in Fig. 3) 

.(0,a/3),, >/\ k a S\y (fe-, k + ) + kpS^y (fc_, fe+) 

JV AA' (K,K)—V F „ • 

(E3) 

Finally, A^y\k,ie, k' ,ie') is the dressed vertex function 
(represented as a triangle in Figs. 3-4), which satisfies 
the self-consistent Bethe-Salpeter equation of Fig. 3c). 
The choice of the quasiparticle self-energy, and the re¬ 
quirement of fulfilling the Ward identities, uniquely de¬ 
termine the form of the Bethe-Salpeter equation, i.e. the 
irreducible interaction I. Fig. 4 shows the three contri¬ 
butions to the irreducible interactions. In formulas, the 
Bethe-Salpeter equation is 

A ( y^\k + ,ie n + ium,k-,ie n ) = A^“ /3) (fc+, fc_) 

3 

+ ^ A ^;^\k + ,ie n + iuj m ,k-,i£ n ) , (E4) 

i—1 

where the three {A y’^\k + ,is n + iu m ,k_,is n ),i = 
1,..., 3} correspond to the three diagrams in fig. 4, and 
read 

A y\ P) (k + ,is n + iu m ,k^,i£ n ) = k B T ^ 

k' ,E n / 

X \\t^iji •> is n ' i^n)G is n ' H - m ) 

+iu m ,k'_,ie n ')G fl (k'_,ie n ') , (E5) 

and 

Ayx P \k +> i£ n + iuJm,k-,i£ n ) = k B T ^ 

k' i£ n ' 

X njjt (k ? ^5 iSri)G/j,' nt 

xA^\k' + ,ie n i +iLJm,k'_ > i£ n ')G ll (k'_,i£ n >) , (E6) 
and finally 

A (k+,ie n + ium,k-,ie n ) = k B T ^ 

k' i £ n' 

x \\t-)k")is n i H - isn + iuJrn)G (fe _|_5 iSfi' + i^m) 
xA^P\k' + ,i£ n '+iUm,k'_,i£„>)G IJ ,(k , _,i£ n >) . (E7) 

Here 

= W(fc - k',iuj m )Vy^(k + ,k' + ) 

xZ> M A(fel,fc_) , (E8) 

where IF(q, iw TO ) is the screened interaction, represented 
as a wavy line in Figs. 3-4. In the large-./Vf limit this is 
given by 




V q 

1 Vq\ nn {Qj IH m ) 


(E9) 


where Xnn(q,u) is the proper density-density response 
function 4 of graphene, ie. 22 

Xnn (q,*w m ) = Nfk B T E E Gx"{q',i £ n ) 

q' ,£n X" ,/i" 

x G M //(q' + q,*£„ + iw m ) 

x iV/HqW + q)^"A"(q / + 9,9') ■ 

(E10) 


Moreover, 

W AA' W '( fc, ’M£n' -*£n) = ^ E 

q',u m f A ",fi" 

xW(q' - q,iu>m' - *w m )£> A 'A"( fc -H fc + - <?') 
xX>A"A(fe+ - <?', k-)D lllt „(k'_ : k' + - q') 
x^'>'(fc) - q',k' + )Gx"{k+ - q',ie n + - iuw) 

xG^"(A;^_ q , i£ n ' -t- ) , (Ell) 

and 

W \\'w'( k '’ + * w m) = ^2 E 

q' X" ,fl" 

xW(q\iuj m f)W(q l - q,iu m > - iuj m )Dxx" {k- , fe_ + q') 
xT>x"X’{k- + q',k + )'D /J ,^(k'_ ) k' + - q') 
x2Vv( fe + - q',k' + )Gx"(k- + q',i£ n + ivm') 
x G [i" (Ay q , i£ n ' T iujm iuj m /) . (E12) 

In Eqs. (E5)-(E12) we suppressed the q-dependence of 
W( 2 \k' ,k,i£ n ’ — i£ n ) and W^ 3 \k',k,i£ n '+ i£ n + iu} m ). 

In what follows we first analytically continue Eqs. (E2)- 
(E12), and then we consider the limit of q = 0, small fre¬ 
quency and low temperature. These limits allows us to 
exactly solve the Bethe-Salpeter equation. This informa¬ 
tion is then used to calculate the low-frequency viscosity 
as shown by Eq. (El). 


1. Analytical continuation to real frequencies 

The analytical continuation of equations similar to 
Eqs. (E2)-(E12) was performed in Ref. 22. There¬ 
fore, we summarize here only the main results. After 
the analytical continuation, the integral on the right- 
hand side of Eq. (E2) contains products of two ad¬ 
vanced Green’s functions (schematically G^G^), two 
retarded ones (G^G^), or one advanced and one re¬ 
tarded (G^G^). The first two combinations have both 
poles on the same side of the complex plane. Therefore, 
in the limit £FT ee 1 they give a negligible contribution, 
and we can retain only the combination G^G^ R \ We 
thus get 

/ d £ 

^[?r F (£ + w) -n F (£)] 

«, A , A ' 

x g[ A) (fc_ , £)A^ (fc_, k + )G^ (fc+, £ + w) 
xA vA ) ( fe +,£++ w > fe —,£-) > 

(E13) 
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where e± = £ ± iO + , and G^\k,e) = G\(k,e+) 

[G^\k,e) = G\ (fe, £—)] is the retarded [advanced] 
Green’s function. The vertex function satisfies the fol¬ 
lowing Bethe-Salpeter equation 


A vf ) ( fe +> e + + w, = A^ ap; (fc+,fc_) 

3 

+ E A va P) (k+,e++w,k-,£-) , 

i =1 


(0 ,a/3), 


(E14) 


where 


A (l,a/3) 

7 i A'A 


de' 


(fc+, £+ + W ,fe_, £ _) = - E / — 


k' 


x [np( £ ') + nB(V - e)] 

x - W&E*'(*'>*><- - £ )' 

x G^ ) (fc+,e' + u) A ^ ) (fe' + ,e' + + u,fe'_,e'_) 

(E15) 


and 


de' 


A^(fc +> e ++W> fc_, e _) = - E I Sn 


k' ,/^,/x' 


x [riF^O + nB(e' - e)] 


Mw (k', k, e'_ -e)- W^, (fe', fe, e' + - e) 


( 2 ) 


X ^V \k+,e' + W ) A y P ^{k' + ,e' + + k'_,s'_) 

xG' A) (fc',e') , 


. («/3) 


and finally 


A 1 


(3,a/3) 
A'A 


de' 


(fc+, £+ + W ,fc_, £ _) = - E / 


k'' 


x [riF(V) + iib(£' + e)] 

x ^IaE(*', k > £ - + £ ) - EaE ( fc '> fe, 4 + £ )' 

xG' A) (fc',£') ■ 


xG™(feV, £ ' + w )A<E(fe' + . £+ + w, feL, O 


(E17) 


In these equations we defined the Fermi and Bose distri¬ 
butions n F (e) = [e E A feBT ) + l] -1 and n B (e) = [e e /( fe BG — 
l] -1 , and the potentials 


( £ - - £ ) - W&w (4 - £ ) = - 4 E 

<j',A ",/x" 

x|W(fe- fe', £ ' — £ )| 2 y [?3 f(ca/ + £ ') — n-p(u/ + e)] 

xSm G^y} (q' — fc,u/ + £ ) 3m G^)(q' — k',u)' + e') 

x Z> A >' (fc + , fe+)U^A(fe-> fe-) 

xlV/Xq' — k,q' — fc')^>/i"A"(g' - feW - fe) , 


(E16) and 


(E18) 


^AAV/i'^- £ ) ^AA'/i/i' ( £ + e ) = E / W+)IE(</,u/_ — w) [?3 f(w' — £ — Uj) — n F (u}' — s' — w)] 


G ( y){k + - q',e + uj — <jo') - G ( y,\k+ — q',e + ui — uj') G^/ 1 (fe+ — q',e' + u — uj') - G[ t t/ ; (fc , + — <?', £ ' + w — uj') 


AA) / 


(RW 


;(A)rfc/ 


xT>A'A"(fe+i fc + - q')'D\"\(k+ - q’, fc_)X> Wi //(fc , _, fc' + - q')'Du.»n>(k' + - q ', fe' + ) , 
and finally 


(E19) 


it; 


(3) 


AA'/x/r 
(R) 


,( £ '_ + £ + w) - w{ a E'( £ + + e + w ) = - E / -w)[n F (w'+ £ )-n F (w'- £ ')] 

/V n" ^ 


G^ ; (fc_ + q', £ + w') - G < 'y){k _ + q', £ + w') G^/(fe+ - q', £ ' + w - a/) - G^f/(fe' + - <?', £ ' + w - <*/) 

xT>aa"( fe-,fe_ + q')'D\"X' (fe- + g',fc+)£> W x"(fe , _,fc+ - q')V ^ (k' + - q',k ' + ) . 


q' ,X" ,/i' 
(R) 


(A), 


(E20) 


2. The Bethe-Salpeter equation 


Here we used that 


Setting q = 0 and taking the limit cu -> 0, Eq. (E13) g (A) (fe) e)G (R) (fe) e + w) ^ 3m[G<? (fe, e)] . 


becomes 


Xctp,[i.v(q — 0 , w ) 


w + j/r e( 


2«w ^-v r de f dn-p{e)\ 

u + i/r ee 2m V 5 £ J 


(E22) 


x3m [G^ r) (fe, £ )] A^ a,3) (fe, fc)EE (fe, £+ + w, fe, £ _) . 

(E21) 
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In writing Eq. (E22) we retained only the singular part 
of the product of the two Green’s functions, i.e. the 
quasiparticle pole, and we neglected the regular part. 
Herein relies our Fermi-liquid approximation. In so do¬ 
ing we are able to significantly simplify the following 
expressions, especially that of the Bethe-Salpeter equa¬ 
tion, which can then be solved analytically. This al¬ 
lows us to determine the viscosity transport time. How¬ 
ever, by neglecting the regular part of the product (E22), 
we miss the Fermi-liquid renormalization of the shear 
modulus Soo- Note that, while the transport time is 
a non-perturbative quantity which diverges in the limit 
a ee —> 0, the renormalization of can be calculated 
perturbatively. Therefore our approach captures the 
most significant effect of e-e interactions in the regime 


in which they are not too strong (i.e. a ee < 1). 

We observe that, at low temperature, the derivative 
of the Fermi function is strongly peaked at e ~ 0. This 
in turn implies that we can set £ = 0 in 3m[G^ R;) (fc,£)] 
on the right-hand side of Eq. (E21). The latter is then 
strongly peaked at the Fermi surface, and allows us to 
set k = k F and A = + in the rest of the integrand. 
As shown in Ref. 22, we cannot set £ = 0 in the fac¬ 
tor A*^(fc,£ + , fc,£_). This happens because the lat¬ 
ter contains Fermi and Bose distributions which combine 
with the factor —dnp(£)/(de) on the right-hand side of 
Eq. (E21) to produce the correct vertex correction,. 

After the analytical continuation Eqs. (E15)-(E20) 
read 


A aa “ /J) (fe,£++w,fc,£_) = - 8 " 

U) + l/T e 


H j Yi / + M e ' -e)][n F (u/ + s') - n F (u' + e)] 

k',q'^ m ** m 


XI W(k- k', s') | 2 3m [Gf (fe', e')] 3m [< G(q' - k, w')] 3m [g™ (q' - fe', u' + £')] V x M (fc, k^V^k', k) 
xT>a >"(<?' — k,q' — k')V^, X n(q' - k',q' - ] (fe', e\ + w, fe', e!_) , 


and 


A * 1 

(2,a/3) 


(E23) 


A AA a/3) (k, £+ + a;, fc, £_) = - 8 ! I Yd I [M e 0 + n B( e ' - e)][n F (w' - e) - n F (u' - e')] 


ds' f dui' 


+ 2mJ 2m 
x \ W(q',u]')\ 2, 3m.[G^(k 1 , e')] 3m[G^(fc — q\ —a/)] 3 m\G^) (k 1 — q 1 , e' — a/) V\\ //(fe, k — q') 
xV y>x (k - q', k)V^(k\ fe' - - q', k')A.^\k',e' + + u>, fe', e'_) , 

and finally 


(E24) 


8 i 


ds 1 f dej' 


A aa ] (k,e + +uj,k,£-) = I Yd I YYi ( n F( £/ )+ «B(e' + £)][n F (w' + e)-n F (u' - e')] 

' k',q' /j,,X" ,fi" 

x |W(g',o; , )| 2 3m[G^(fe , ,£ , )]3m[G^(fe + q',w')]3m[G^(fe' - q',s' - uj')]T>x\"(k, k + q') 


xV x „ x (k + q',k)V^„(k',k' - g')2V>(fc' - q',k')A<fif\k',e' + + u,k' ,e'_) . 


(E25) 


In these equations the limit to —► 0 is understood. 


In evaluating an integral of the form 


We now recall that Eqs. (E23)-(E25) must be intro¬ 
duced into Eq. (E21) and integrated over e with the 
weighting factor dn F (£)/ (ds). To perform such integra¬ 
tion we use the fact that 


. dlZ F (E^ r / /\ / / ,1 

J\f = —[n F (s ) + n B (£ - e)] 
x [ n F (s" + e) — n F (£" + e')] 

= de „ (e + E ) -« F (e)] 

x \n F (e' + e") — n F (£ r )] 
„ 2 driB(E") dn F (£) dn F (E r ) 
ds" ds ds' 


(E26) 


1 = r d£"^P-£" 2 f(£") , (E27) 

where f(s") is some smooth function of its argu¬ 
ment, we exploit the fact that the weighting function 
£" 2 <9n B (£")/<9£" is strongly peaked at e" = 0 and its 
width scales with fc|T 2 /E F . This does not mean, how¬ 
ever, that one can simply replace f(e") by /(0). Such 
a crude approximation may introduce spurious diver¬ 
gences 22 . To take this into account we approximate 

J= _ 27r-(fc B r)- / (-) +0(T 4 ) , ( E28 ) 

where 22 £ = ^fc B T is approximated with half the variance 
of the distribution e 2 dn^s) / (ds). Therefore, ( = ir/x/b. 
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With this approximation we can rewrite 

4 i{k B T) 2 


A 


AA 


(fc,u/ + , fc, 0 ) = — 


E E 


3(w + */r ee ) „ 

v ' k' ,q' n,fi" 


X| W(fc - fc', 0) | 2 9m [G^ (fc', 0)] Sm [g$ } ( g' - fc, 0) 


xSm 


G^V-fc'.O) V Xfl (k,k')V^(k',k) 


xI?A'>"(g' — k,q' — k'yD^’x”^ - fc',g' - fc) 
xA^(fc> + ,fc',(r) , (E29) 


and 


A 1 


(2, a/3) 
AA 


(fc,w+,fc,0 )= - 


4i(fc B T) 2 

3(w + t/ree)^ /;iiA „ iM „ 


E E 


x|W(g',0)| 2 3m[G( t R) (fc , ,0)]3m[G^ i : ) (fc-g , ,0)] 

xSm[Gg } (k' - g', 0)] P AA »(fc, fc - q')V x „ x (k - g', k) 

xX> MA1 "(fc',fc' - gO^Vv^ - <?'. fc ') 
xA^(k',u} + ,k',0-) , (E30) 

and finally 

4i(fc B T) 2 


A^\k, U+ ,k,0 ) = 


3 (w + f/'T’ee) 


E E 


k\q' fi,X" ,fi" 

x|W(g',0)| 2 STO[G( t R) (fc',0)]Sm[G^ ) (fc + g',0)] 
x3m[GW(fe' - g', 0)]Z> AA »(fc, k + q')V x „ x {k + g', fc) 
xH w "(fc', k' - q')V^n^(k' - g', fc') 
xA^f (fc',oc + ,fc',(T) . (E31) 

Shifting fc' —>- fc — g" and g' —>■ fc — fc" in Eq. (E30), we 
obtain 


A 


4i(fc B T) 2 


£ a/s) (fc,u; + ,fc,0 ) = - , —r 

3(W + l / T e e) 


E E 

k",q" fj.,X" ,/j." 


x|W(fc-fc",0)| 2 3TO[G( 1 R) (fc-g",0)]3m[G^ ) (fc",0)] 
x9m[G<*>(fc" - g", 0)1 V xx „ (fc, fc")2V A (fc", fc) 


xV(fc - g",fc" - g")2W*" - 9",fe - g") 

xA^f } (fc -g",w + ,fc-g",0 _ ) , 


(E32) 


while shifting fc' —► fc" —g" and g' —)• fc" —fc in Eq. (E31), 
we get 


A£«»0fe,„ + ,M-)= 4i( ' iBT) 


E E 


3(w + */Tee) „ 

x| W{k" - fc,0)| 2 9?n[G( t R) (fc" - g", 0)]3m[G™(fc", 0)] 
x9m[G{ t R) (fc - g", 0)]Z> AA »(fc, fe")2V A (fe", fc) 

xV(fc" - - g'OzW* - g",*" - g") 

xA(“«(fc"-g",a; + ,fc"-g",0-) . (E33) 

We now observe that the products of the three Green’s 
functions on the right-hand sides of Eqs. (E29)-(E31) 


constrains (i) /r = n" = A" = + and (ii) the momenta 
of their arguments to be at the Fermi surface. After 
renaming dummy momentum variables and noting that 
£>aa '(-fc, -k’) = V xy (k, fc'), we finally get 

X E E |W / (fc-fc',0)| 2 9m[GW(fc',0)] 

k' ,q' ,X' f 


G ( y) (g' - fc, 0)J 9m [G^ j (g' - fc', 0) 
xV Xfl (k, fc')2? MA (fc', fc)2A A " M "(g' - fc, g' - fc') 
x£V" A "(g' - fc',g' - fc)[A^f ) (fc',w + ,fc',0“) 
+A^f } (fc - g',w+,fc - g',0~) 

—A^f } (fc' - g',o; + ,fc' - g',0 - )] . 

We solve Eq. (E34) with the following Ansatz: 


-.(H)/ 


(E34) 


A^ ) (fc,cj + ,fc,0 ) = 7(w)A RJ i | a/3J (fc,fc) , (E35) 

which reduces the self-consistent equation to an algebraic 
one. We then project the result along A++^\k,k), i.e. 

we multiply it by A^°(f 7 ^(fc, fc) and we take the trace over 
the greek indices. We get 

4i(fc B T) 2 


(0,a/3), 


y(w) = 1 - 7(w) 


3(w + i/r ee ) “' 

Ac ,q /z,/r ,A 


E E ime(/c - fc',o)i= 


x9m[G^ R) (fc , ,0)]9mfG^ ) (g' - fc, 0) 


(E36) 


G^V - fc', 0)] 2? v (fc, fc')^ A (fc', fc) 
x£> A " M "(g' -k,q' - fc')X> M -/ A //(g' - fc',g' - fc) 

X [cOS 2 (</?fc/) + COS 2 (v2fc_g.) - C0S 2 (v3fc-_q/)] , 

Here we used that k = k' = kp and that the Green’s func¬ 
tions of Eq. (E34) constrain their momentum arguments 
to be at the Fermi surface. We also assumed fc = kx. 
Eq. (E36) can be rewritten as 

i/ r ee - i/r v 


7M = 1 +7(w)- 


+ i / 7 


(E37) 


where 

2. = 4(fc B T) 2 
t v 3 


E E \w(k-k\o)f 

k' ,q' ,X" 


x9m[G{ 1 R) (fc',0)]9m[G®(g' - fc,0) 

x9m [g< r) ( g' - fc', 0)] V X/1 (k, fc')^ A (fc', fc) 
x!> A " M "(g' fc,g' - k')D fJ ," X "( g' - fc',g' - fc) 

X [cOS 2 (^fc.) + COS 2 (^fc„q') - COS 2 (v3fc'_q') - l] . 

(E38) 

Eq. (E37) is readily solved by 

w + f/r ee 


7(w) = 


oj + i/r v 


(E39) 
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which contains all the information about the vertex cor¬ 
rections. 

Finally, putting Eqs. (E35) and (E39) back into 
Eq. (E13) and Eq. (El), and taking the limit w —> 0, 
we get 


derived in a regime in which intraband transitions are 
responsible for the dominant contribution, in what fol¬ 
lows we approximate the formulas of Ref. 5 by neglecting 
contributions due to interband processes. This approxi¬ 
mation is valid in the low-temperature regime. We get 


* v 

= — E 

ee fc,A,A' 


de dn F (e) ( 


2-7T 

xG^(k,s) a[°; x 


de 

\k,k) 


G[ A \k,e)A. 


a 7 ] (k,k) 




Y Sm [g^ A) {k, 0)G™ (fc, 0)1 [a^ (fc, fc) 1 2 


1 

= -ne F r v . 
4 


(E40) 


One immediately recognizes So = ne p/4 as the non¬ 
interacting bulk modulus (which is not renormalized by 
e-e interactions, as discussed before) and r v as the vis¬ 
cosity transport time. 


1 4 [°° dn F (f) [ + °° 1 -n F «-o;) 

r v - (2 tt) 2 J * dt; J W 1 - exp(—/3w) 


n -)-oo 


dq q 


e{q,u,T) 

1 - 


x A ++ (k F ,q,u) 

x 4 (q 2 -u 2 /v 2 ) ^ 
2fcp(fcp — <jj/v f) 


^rn[x { n^(q,uJ,T)] 


q 2 — u 2 /i 


4k F (k F — oj/v f) 


(E43) 


where 


3. Transformation of Eq. (E38) to a 
computationally efficient formula 


We start from the definition of Eq. (E38), and we shift 
k! -» fc — q, q' —> k" -f fc. We get 

- = ^P^Y E \W( qi 0)\ 2 Zm[GW(k-q,0)} 


G™(fc",0) 


G™(k" + q, 0) 


x T>\ M (fc, fc - q)V^\(k - q, k)T>x»^'{k", k" + q) 

X 22/j"A" (fc” + 9: fc”) 

X [cOS 2 ((/?fc_g) +COS 2 (</5fc//) - C0S 2 (v5fc// + q) - l] . 

(E41) 


Now we use the fact that, if the scattering occurs at the 
Fermi surface, k" + q (fe") is opposite to fc (fc + q) (see 
also Ref. 22). Therefore 


1 

Tv 


8{k B T) 2 


E E 

q.q' 


|W(q,0)| 2 3m G™(fc",0) 


x 3m [Gf) (fc — q, 0)] 3m 


G^(fc" 


9,0) 


x V x M (fc, fc - q)D^\(k - q, fc)X>A'>"(fc", fc" + q) 
x V^x"(k" + q. fc")[l - cos 2 (^ fc+ q)] . (E42) 


Note that the expression of Eq. (E42), if one exclude the 
matrix element in its last line, coincides with the quasi¬ 
particle lifetime for k B T <C £f reported in Refs. 21 and 
22. The expression of the quasiparticle lifetime valid at 
all temperatures was also given in Ref. 5. At low tem¬ 
perature (a limit that almost always holds in graphene), 
we can thus rewrite l/r v using the formulas provided 
in Ref. 5, amending them by multiplying the integrand 
with the matrix element 1 — cos 2 (ipk +q ). Since that was 


4(fc — ui/v f) 

v F x/\(2k - u/v F ) 2 - q 2 ](q 2 - w 2 /^f) 
x ^ g 2 -co 2 /v 2 - 

4 fc(fc — u/v f) 

x ©| [(2fc - u/v F ) 2 - q 2 ] (q 2 - w 2 /u|)} , 

(E44) 


and we used that 


cos ((pfc+g) y 1 


q 2 - co 2 /v 2 
2 k(k — u/v f) 


(E45) 


Eq. (E43) has a form similar to the quasiparticle con¬ 
tribution to the decay rate given in Ref. 5. We did 
not write the quasihole contribution, but we included it 
with an extra factor of two that multiplies the whole 
expression. Indeed, in the low-temperature limit the 
quasiparticle and quasihole contributions are identical 5 . 
The integration over e in Eq. (E43), with the weight¬ 
ing factor <9n F (£)/(<9£), is the same integration we per¬ 
formed in the Bethe-Salpeter equation [see the discus¬ 
sion after Eq. (E25)]. Its origin has to be found in the 
fact that, in the limit of low temperature, the energy- 
dependent self-energy can be replaced by its average over 
the thermally excited states. Note that, since in the low- 
temperature limit the function dn F (£)/(c?£) is strongly 
peaked at £ = 0 (which in turn implies k = k F ), we set 
k = k F everywhere in the integrand on the right-hand 
side of Eq. (E43) , except in the Fermi function on its first 
line, which strongly depends on £. 

Following Refs. 5, 21, and 22, we can now calculate 
Eq. (E43). As noted in the main text, the matrix element 
1 — cos 2 (<pfc +(J ) completely suppresses the divergence due 
to the small-q region. Following the same manipulations 
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performed in Ref. 5 we get 


straightforward and the result reads (restoring H) 


- = 327V(0)a 2 e [ [ dto un B (-«) 

Tv J — oo J —oo 

«f(£ - w) 


/»2/Cp —CD / 


G?g 


■'MA’f Q 

1 ~ g 2 /(4fcg) 

g 2 — oj 2 /vp 
q 2 - Ul 2 /v 2 
2k(kp — uj/vp) 


1 + 


1 - 


9tf ^ g| F ^ 2 1 - g 2 /(4fcg) 

g 2 Up g 2 — uj 2 /v p 
q 2 - U 2 /v 2 p 


4fcp(fcp — ui/vp) 


(E46) 


Here N( 0) = 7Vf£F/[27r(?i'UF) 2 ] is the density-of-states at 
the Fermi energy. Since the ratio in the integrand on the 
first line of Eq. (E46) is strongly peaked around £, w = 0, 
we take the limit of £, u> —> 0 in the rest of the inte¬ 
grand. Contrary to the case of the quasiparticle lifetime 
this does not lead to any divergence. Indeed, the inte¬ 
grand is regular for q —1 0, 2kp. The integration is then 


1 = -^%41! ai [42 - 52JV,„« - 15A '?& 
r v 9 he p 

+ 157V f 3 a 3 e - 3(16 - 247V 2 a 2 e + 5 Nfa^) 

X COth~ 1 (l + A^fOee)] ■ 

(E47) 

To get this equation we used that gxF = NfOt ee kp, and 
that in the low temperature limit 


dt ; 


dnp(£,) 


di . 

In the limit of a, 


[ +co 1 -n F (g-a;) 
I -oo 1 - exp(-/3w) 


0 we get 


37T 2 

16 


(k B T) 2 . 

(E48) 


1 

T v 


-> -N { 


87T (fc B T) 2 
3 hep 


a 2 e ln(a ee ) • 


(E49) 
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